空汙資料庫期末報告(G03)
作者:楊淳光 Oscar
2023-05-12 09:57
1.資料介紹與下載
1.1資料來源:
1.2資料下載:
-
利用清單關鍵字收尋：「空氣品質監測小時值(一般污染物,每日更新)」
選擇代碼為AQX_P_13
下滑至歷史資料下載區塊
勾選要下載的年-月檔案
點選「下載勾選項目」，選擇csv檔(格式)
下載後解壓縮得到csv 格式的檔案
利用以下程式碼可將檔案匯入RStudio裡
- 設定工作區域
setwd("E:/rumiworkspace")- 讀取資料
filename = "./data/G03_P13_2014_2022/空氣品質監測小時值資料(一般污染物,每日更新) (2014-01).csv"
df <- read.csv(filename,fileEncoding = "UTF-8-BOM")- 檢查資料
head(df)## X.siteid. X.sitename. X.itemid. X.itemname. X.itemengname. X.itemunit.
## 1 11 士林 5 氮氧化物 NOx ppb
## 2 11 士林 6 一氧化氮 NO ppb
## 3 25 頭份 33 細懸浮微粒 PM2.5 μg/m3
## 4 25 頭份 14 溫度 AMB_TEMP ℃
## 5 11 士林 7 二氧化氮 NO2 ppb
## 6 3 萬里 11 風向 WIND_DIREC degrees
## X.monitordate. X.monitorvalue00. X.monitorvalue01. X.monitorvalue02.
## 1 2014-01-01 16 15 29
## 2 2014-01-01 1.2 .8 3.3
## 3 2014-01-01 57 62 51
## 4 2014-01-01 14 13 13
## 5 2014-01-01 15 14 25
## 6 2014-01-01 212 210 213
## X.monitorvalue03. X.monitorvalue04. X.monitorvalue05. X.monitorvalue06.
## 1 86 89 93 101
## 2 43 49 56 64
## 3 40 33 33 33
## 4 12 12 12 11
## 5 43 40 37 36
## 6 212 213 213 211
## X.monitorvalue07. X.monitorvalue08. X.monitorvalue09. X.monitorvalue10.
## 1 119 129 20.0 14.0
## 2 82 82 5.7 3.3
## 3 38 39 34.0 31.0
## 4 11 14 17.0 20.0
## 5 37 47 15.0 11.0
## 6 209 192 69.0 44.0
## X.monitorvalue11. X.monitorvalue12. X.monitorvalue13. X.monitorvalue14.
## 1 11 16 12 15
## 2 2.9 3.8 2.9 2.7
## 3 30 35 37 35
## 4 21 22 21 22
## 5 8.5 12 9.5 12
## 6 58 49 50 48
## X.monitorvalue15. X.monitorvalue16. X.monitorvalue17. X.monitorvalue18.
## 1 16 14 17 14
## 2 2.5 2.3 1.9 1.9
## 3 38 45 51 59
## 4 22 21 20 19
## 5 14 12 15 12
## 6 57 77 66 200
## X.monitorvalue19. X.monitorvalue20. X.monitorvalue21. X.monitorvalue22.
## 1 13 11 11 8.6
## 2 1.7 1.9 1.9 1.6
## 3 65 62 56 51.0
## 4 18 17 17 17.0
## 5 11 8.7 9 6.9
## 6 208 205 206 200.0
## X.monitorvalue23.
## 1 6.8
## 2 1.7
## 3 47.0
## 4 17.0
## 5 5.1
## 6 207.01.3變項介紹:
•SiteId(測站代碼):
每個測站都有一個自己的代碼，代表該測站的識別(1~72、75、77、78、80、83、84、85)
•SiteName(測站名稱)：
是測站的名稱(基隆、汐止、萬里、…、大城)
•ItemId(測項代碼)：
每個測項都有一個自己的代碼，代表該測項的識別(1~9、11、14、21~23、31~33、36、38、143~144)
•ItemName(測項名稱)：
是測項的名稱(二氧化硫、一氧化碳、臭氧、…、小時風向值)
•ItemEngName(測項英文名稱)：
是測項的英文名稱(SO2、CO、O3、…、WD_HR)
•ItemUnit(測項單位)：
測項的計量單位，如 μg/m³、ppm等。
Item表格整理+說明✔
| ItemId | ItemName | ItemEngName | ItemUnit | 測項說明 |
|---|---|---|---|---|
| 1 | 二氧化硫 | SO2 | ppb | 空氣中的一種污染物，主要來自燃煤、石油等化石燃料的燃燒排放。 |
| 2 | 一氧化碳 | CO | ppm | 空氣中的一種有毒氣體，主要來自車輛排放和燃燒過程中不完全燃燒。 |
| 3 | 臭氧 | O3 | ppb | 空氣中的一種臭氣，由氮氧化物和揮發性有機物在陽光照射下產生，對人體和植物有害。 |
| 4 | 懸浮微粒 | PM10 | μg/m3 | 空氣中的固體和液體顆粒物，可造成呼吸道、心血管等健康問題。 |
| 5 | 氮氧化物 | NOx | ppb | 包括一氧化氮和二氧化氮，主要來自燃煤、汽車等排放源。 |
| 6 | 一氧化氮 | NO | ppb | 一種有毒氣體，主要來自汽車排放和燃燒過程中高溫氮氧化物的生成。 |
| 7 | 二氧化氮 | NO2 | ppb | 一種有毒氣體，主要來自汽車、工廠等排放源，可造成呼吸道問題。 |
| 8 | 總碳氫化物 | THC | ppm | 包括甲烷、乙烷、丙烷等氣體，以及揮發性有機化合物，主要來自燃料燃燒和工業過程。 |
| 9 | 非甲烷碳氫化合物 | NMHC | ppm | 不包括甲烷的碳氫化合物，主要來自石油和天然氣開采、運輸等過程。 |
| 11 | 風向 | WIND_DIREC | degrees | 風的方向。 |
| 14 | 溫度 | AMB_TEMP | ℃ | 空氣的溫度，通常用攝氏度或華氏度表示。 |
| 21 | 酸雨 | PH_RAIN | pH | 指大氣中酸性物質與雨水混合而成的雨水，對生態環境和建築物有害。 |
| 22 | 導電度 | RAIN_COND | μmho/cm | 水體中離子的濃度和移動能力的度量，通常用西門子/厘米表示。 |
| 23 | 雨量 | RAINFALL | mm | 一定時間內的雨水總量，通常用毫米表示。 |
| 31 | 甲烷 | CH4 | ppm | 一種無色、無味、易燃的氣體，是溫室氣體之一。 |
| 32 | 降雨強度 | RAIN_INT | ㎜ | 降雨時單位時間內雨量的多寡，常用毫米每小時（mm/hr）來表示。 |
| 33 | 細懸浮微粒 | PM2.5 | μg/m3 | 空氣中直徑小於2.5微米的顆粒物，會對人體健康造成危害，如引發呼吸道疾病。 |
| 36 | 二氧化碳 | CO2 | ppm | 一種無色、無味、不易燃的氣體，是主要的溫室氣體之一，對全球氣候變化有重要影響。 |
| 38 | 相對濕度 | RH | % | 空氣中水蒸氣的含量與空氣中最大可能的含量之比，通常以百分比表示。 |
| 143 | 小時風速值 | WS_HR | m/sec | 一小時內風速的最大值。 |
| 144 | 小時風向值 | WD_HR | degrees | 一小時內風向的平均值。 |
•MonitorDate(監測日期)：
測量數據的日期(20xx/xx/xx)
•MonitorValue0(0時數值)~MonitorValue23(23時數值)：
測量數據的具體數值，代表在該時的空氣品質監測值，單位為測項的計量
1.4變項之間相關性介紹:
SiteId 和 SiteName 是相關的，表示同一個測站。
在分析數據時，可以根據 SiteName 篩選出特定城市或地區的監測數據，或者根據 SiteId 查詢該測站的歷史數據。
MonitorDate 和 MonitorValue0~MonitorValue23 是相關的，表示同一天的不同時間點的監測數據。
可以分析數據的時間趨勢和變化情況。
ItemId 和 ItemName/ItemEngName：每個 ItemId 都對應一個 ItemName 和一個 ItemEngName，表示測量的物質名稱和英文名稱，因此這三個變項是相關的。
example✔
在分析數據時，可以根據 ItemName/ItemEngName 篩選出特定物質的監測數據，或者根據 ItemId 查詢該測項的歷史數據。
setwd("E:/rumiworkspace")
filename = "./data/datatest/空氣品質監測小時值資料(一般污染物,每日更新)_test.csv"
df = read.csv(filename,fileEncoding = "UTF-8-BOM")
datatable(head(df, 100), options = list(scrollx = TRUE, scrollCollapse = TRUE))2.資料整理
目標 利用 R 的 dplyr 的套件中的指令，來整理資料，並且轉置資料。
2.1 了解並載入dplyr套件
#載入dplyr套件
library(dplyr)
補充：淺談dplyr基礎功能
select()：從數據中選擇列
filter()：數據行的子集 (篩選符合條件的觀察值)
group_by()：匯總數據 (根據資料內的類別變數對另一個變數進行統計)
summarise()：匯總數據(計算匯總統計信息)
arrange()：排序數據 (變數遞增或是遞減排序資料)
mutate()：創建新變量
示範：
#自訂新變量[newcode]
setwd("E:/rumiworkspace")
airpollution<-read.csv('./data/datatest/空氣品質監測小時值資料(一般污染物,每日更新)_test.csv'
,fileEncoding = "UTF-8-BOM")
airpollution.part<-mutate(airpollution, newcode=round(as.numeric(averageday)))
datatable(head(airpollution.part, 100), options = list(scrollx = TRUE, scrollCollapse =
TRUE))#利用filter可以抓取指定測項和其規範數值
setwd("E:/rumiworkspace")
airpollution<-read.csv('./data/datatest/空氣品質監測小時值資料(一般污染物,每日更新)_test.csv'
,fileEncoding = "UTF-8-BOM")
airpollution.filter<-filter(airpollution, itemengname == 'NOx',as.numeric(averageday) >
100)
airpollution.filter## itemname itemengname itemunit monitordate averageday
## 1 氮氧化物 NOx ppb 2014/1/1 135.63
## 2 氮氧化物 NOx ppb 2014/1/2 110.00
## 3 氮氧化物 NOx ppb 2014/1/2 194.73
## 4 氮氧化物 NOx ppb 2014/1/2 100.98
## 5 氮氧化物 NOx ppb 2014/1/2 129.25
## 6 氮氧化物 NOx ppb 2014/1/3 174.67
## 7 氮氧化物 NOx ppb 2014/1/3 106.13
## 8 氮氧化物 NOx ppb 2014/1/4 142.96
## 9 氮氧化物 NOx ppb 2014/1/5 113.38
## 10 氮氧化物 NOx ppb 2014/1/6 124.17
## 11 氮氧化物 NOx ppb 2014/1/7 130.46
## 12 氮氧化物 NOx ppb 2014/1/7 102.29
## 13 氮氧化物 NOx ppb 2014/1/7 193.88
## 14 氮氧化物 NOx ppb 2014/1/8 110.54
## 15 氮氧化物 NOx ppb 2014/1/8 187.17
## 16 氮氧化物 NOx ppb 2014/1/9 124.88
## 17 氮氧化物 NOx ppb 2014/1/11 119.13
## 18 氮氧化物 NOx ppb 2014/1/12 124.882.2 資料察看
以了解我們資料是否有地方可以修改。
# 檔案名稱
setwd("E:/rumiworkspace")
filename = "./data/G03_P13_2014_2022/空氣品質監測小時值資料(一般污染物,每日更新) (2014-01).csv"
# 讀入資料
df = read.csv(filename, fileEncoding="UTF-8-BOM")
str(df)## 'data.frame': 27756 obs. of 31 variables:
## $ X.siteid. : int 11 11 25 25 11 3 42 25 71 11 ...
## $ X.sitename. : chr "士林" "士林" "頭份" "頭份" ...
## $ X.itemid. : int 5 6 33 14 7 11 4 11 10 10 ...
## $ X.itemname. : chr "氮氧化物" "一氧化氮" "細懸浮微粒" "溫度" ...
## $ X.itemengname. : chr "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## $ X.itemunit. : chr "ppb" "ppb" "μg/m3" "℃" ...
## $ X.monitordate. : chr "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
## $ X.monitorvalue00.: chr "16" "1.2" "57" "14" ...
## $ X.monitorvalue01.: chr "15" ".8" "62" "13" ...
## $ X.monitorvalue02.: chr "29" "3.3" "51" "13" ...
## $ X.monitorvalue03.: chr "86" "43" "40" "12" ...
## $ X.monitorvalue04.: chr "89" "49" "33" "12" ...
## $ X.monitorvalue05.: chr "93" "56" "33" "12" ...
## $ X.monitorvalue06.: chr "101" "64" "33" "11" ...
## $ X.monitorvalue07.: chr "119" "82" "38" "11" ...
## $ X.monitorvalue08.: chr "129" "82" "39" "14" ...
## $ X.monitorvalue09.: chr "20.0" "5.7" "34.0" "17.0" ...
## $ X.monitorvalue10.: chr "14.0" "3.3" "31.0" "20.0" ...
## $ X.monitorvalue11.: chr "11" "2.9" "30" "21" ...
## $ X.monitorvalue12.: chr "16" "3.8" "35" "22" ...
## $ X.monitorvalue13.: chr "12" "2.9" "37" "21" ...
## $ X.monitorvalue14.: chr "15" "2.7" "35" "22" ...
## $ X.monitorvalue15.: chr "16" "2.5" "38" "22" ...
## $ X.monitorvalue16.: chr "14" "2.3" "45" "21" ...
## $ X.monitorvalue17.: chr "17" "1.9" "51" "20" ...
## $ X.monitorvalue18.: chr "14" "1.9" "59" "19" ...
## $ X.monitorvalue19.: chr "13" "1.7" "65" "18" ...
## $ X.monitorvalue20.: chr "11" "1.9" "62" "17" ...
## $ X.monitorvalue21.: chr "11" "1.9" "56" "17" ...
## $ X.monitorvalue22.: chr "8.6" "1.6" "51.0" "17.0" ...
## $ X.monitorvalue23.: chr "6.8" "1.7" "47.0" "17.0" ...head(df)## X.siteid. X.sitename. X.itemid. X.itemname. X.itemengname. X.itemunit.
## 1 11 士林 5 氮氧化物 NOx ppb
## 2 11 士林 6 一氧化氮 NO ppb
## 3 25 頭份 33 細懸浮微粒 PM2.5 μg/m3
## 4 25 頭份 14 溫度 AMB_TEMP ℃
## 5 11 士林 7 二氧化氮 NO2 ppb
## 6 3 萬里 11 風向 WIND_DIREC degrees
## X.monitordate. X.monitorvalue00. X.monitorvalue01. X.monitorvalue02.
## 1 2014-01-01 16 15 29
## 2 2014-01-01 1.2 .8 3.3
## 3 2014-01-01 57 62 51
## 4 2014-01-01 14 13 13
## 5 2014-01-01 15 14 25
## 6 2014-01-01 212 210 213
## X.monitorvalue03. X.monitorvalue04. X.monitorvalue05. X.monitorvalue06.
## 1 86 89 93 101
## 2 43 49 56 64
## 3 40 33 33 33
## 4 12 12 12 11
## 5 43 40 37 36
## 6 212 213 213 211
## X.monitorvalue07. X.monitorvalue08. X.monitorvalue09. X.monitorvalue10.
## 1 119 129 20.0 14.0
## 2 82 82 5.7 3.3
## 3 38 39 34.0 31.0
## 4 11 14 17.0 20.0
## 5 37 47 15.0 11.0
## 6 209 192 69.0 44.0
## X.monitorvalue11. X.monitorvalue12. X.monitorvalue13. X.monitorvalue14.
## 1 11 16 12 15
## 2 2.9 3.8 2.9 2.7
## 3 30 35 37 35
## 4 21 22 21 22
## 5 8.5 12 9.5 12
## 6 58 49 50 48
## X.monitorvalue15. X.monitorvalue16. X.monitorvalue17. X.monitorvalue18.
## 1 16 14 17 14
## 2 2.5 2.3 1.9 1.9
## 3 38 45 51 59
## 4 22 21 20 19
## 5 14 12 15 12
## 6 57 77 66 200
## X.monitorvalue19. X.monitorvalue20. X.monitorvalue21. X.monitorvalue22.
## 1 13 11 11 8.6
## 2 1.7 1.9 1.9 1.6
## 3 65 62 56 51.0
## 4 18 17 17 17.0
## 5 11 8.7 9 6.9
## 6 208 205 206 200.0
## X.monitorvalue23.
## 1 6.8
## 2 1.7
## 3 47.0
## 4 17.0
## 5 5.1
## 6 207.0從以上兩表發現有幾處是需要修改的。
- 欄位名稱需修改
| 原名稱 | 英文名稱 | 中文名稱 |
|---|---|---|
| X.siteid. | siteid | 測站代碼 |
| X.sitename. | sitename | 測站名稱(V) |
| X.itemid. | itemid | 測項代碼 |
| X.itemname. | itemname | 測項名稱 |
| X.itemengname. | itemengname | 測項英文名稱(V) |
| X.itemunit. | itemunit | 測項單位 |
| X.monitordate. | monitordate | 測量日期(V) |
| X.monitorvalue00. ~ 23 | h00~h23 | 測量值(V) |
欄位資料類型需要修改
monitordate須為日期資料類型測量值需為數值資料類型
2.3 測站資訊
查看測站相關變數。
有下列資訊：
測站代碼(siteid)
測站名稱(sitename)
# 測站資訊檔
stateInfo = data.frame(df$X.siteid., df$X.sitename.)
colnames(stateInfo) = c("測站代碼", "測站名稱")
stateInfo = unique(stateInfo)
str(stateInfo)## 'data.frame': 76 obs. of 2 variables:
## $ 測站代碼: int 11 25 3 42 71 35 56 72 12 4 ...
## $ 測站名稱: chr "士林" "頭份" "萬里" "嘉義" ...2.4 測項資訊
查看測項相關變數。
有下列資訊：
測項代碼(itemid)
測項中文名稱(itemname)
測項英文名稱(itemengname)
測項單位(itemunit)
# 測項資訊檔
itemInfo = data.frame(df$X.itemid., df$X.itemname.,
df$X.itemengname., df$X.itemunit.)
colnames(itemInfo) = c("測項代碼", "測項中文名稱", "測項英文名稱", "測項單位")
itemInfo = unique(itemInfo)
str(itemInfo)## 'data.frame': 12 obs. of 4 variables:
## $ 測項代碼 : int 5 6 33 14 7 11 4 10 3 38 ...
## $ 測項中文名稱: chr "氮氧化物" "一氧化氮" "細懸浮微粒" "溫度" ...
## $ 測項英文名稱: chr "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## $ 測項單位 : chr "ppb" "ppb" "μg/m3" "℃" ...2.5 測量資料(欄位名稱)更改
以便於我們識別欄位。
# 更改欄位名稱
dataInfo = dplyr::select(
df, X.sitename., X.itemengname.,
X.monitordate., starts_with("X.monitorvalue")
)
colnames(dataInfo) = c("測站名稱", "測項英文名稱", "測量日期",
paste("h0",seq(0,9),sep=""),
paste("h",seq(10,23),sep=""))
str(dataInfo)## 'data.frame': 27756 obs. of 27 variables:
## $ 測站名稱 : chr "士林" "士林" "頭份" "頭份" ...
## $ 測項英文名稱: chr "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## $ 測量日期 : chr "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
## $ h00 : chr "16" "1.2" "57" "14" ...
## $ h01 : chr "15" ".8" "62" "13" ...
## $ h02 : chr "29" "3.3" "51" "13" ...
## $ h03 : chr "86" "43" "40" "12" ...
## $ h04 : chr "89" "49" "33" "12" ...
## $ h05 : chr "93" "56" "33" "12" ...
## $ h06 : chr "101" "64" "33" "11" ...
## $ h07 : chr "119" "82" "38" "11" ...
## $ h08 : chr "129" "82" "39" "14" ...
## $ h09 : chr "20.0" "5.7" "34.0" "17.0" ...
## $ h10 : chr "14.0" "3.3" "31.0" "20.0" ...
## $ h11 : chr "11" "2.9" "30" "21" ...
## $ h12 : chr "16" "3.8" "35" "22" ...
## $ h13 : chr "12" "2.9" "37" "21" ...
## $ h14 : chr "15" "2.7" "35" "22" ...
## $ h15 : chr "16" "2.5" "38" "22" ...
## $ h16 : chr "14" "2.3" "45" "21" ...
## $ h17 : chr "17" "1.9" "51" "20" ...
## $ h18 : chr "14" "1.9" "59" "19" ...
## $ h19 : chr "13" "1.7" "65" "18" ...
## $ h20 : chr "11" "1.9" "62" "17" ...
## $ h21 : chr "11" "1.9" "56" "17" ...
## $ h22 : chr "8.6" "1.6" "51.0" "17.0" ...
## $ h23 : chr "6.8" "1.7" "47.0" "17.0" ...2.6 轉置資料(寬表與長表)
何謂轉置?
在線性代數的矩陣操作裡，依據左上至右下的對角線來進行列資料與欄資料的對換稱為「轉置」（Transpose）；在表格式資料（試算表、關聯表格或資料框）中也有相似的操作，由於表格式資料相對於矩陣，多出了變數的欄位名稱與觀測值的列索引，轉置前後的外觀又被稱作寬表（Wide Format）與長表（Long Format）。
為何要轉置資料呢?
資料的外觀與函式並無法呼應。這時資料轉置的需求就相當地明確
我們有兩種欄位類型：
- 長表(long)：測站名稱，測項名稱，測量日期，測量小時，測量值
# 轉製成長版
varying = c(paste("h0",seq(0,9),sep=""),
paste("h",seq(10,23),sep=""))
df.long = reshape(dataInfo, varying=varying,
v.names="測量值",
timevar="測量小時",
#times=varying, #seq_along(varying[[1]]),
direction="long")
str(df.long)## 'data.frame': 666144 obs. of 6 variables:
## $ 測站名稱 : chr "士林" "士林" "頭份" "頭份" ...
## $ 測項英文名稱: chr "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## $ 測量日期 : chr "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
## $ 測量小時 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ 測量值 : chr "16" "1.2" "57" "14" ...
## $ id : int 1 2 3 4 5 6 7 8 9 10 ...
## - attr(*, "reshapeLong")=List of 4
## ..$ varying:List of 1
## .. ..$ 測量值: chr [1:24] "h00" "h01" "h02" "h03" ...
## .. ..- attr(*, "v.names")= chr "測量值"
## .. ..- attr(*, "times")= int [1:24] 1 2 3 4 5 6 7 8 9 10 ...
## ..$ v.names: chr "測量值"
## ..$ idvar : chr "id"
## ..$ timevar: chr "測量小時"df.long$id = NULL
attr(df.long, "reshapeLong") = NULL
str(df.long)## 'data.frame': 666144 obs. of 5 variables:
## $ 測站名稱 : chr "士林" "士林" "頭份" "頭份" ...
## $ 測項英文名稱: chr "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## $ 測量日期 : chr "2014-01-01" "2014-01-01" "2014-01-01" "2014-01-01" ...
## $ 測量小時 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ 測量值 : chr "16" "1.2" "57" "14" ...head(df.long, n=30)## 測站名稱 測項英文名稱 測量日期 測量小時 測量值
## 1.1 士林 NOx 2014-01-01 1 16
## 2.1 士林 NO 2014-01-01 1 1.2
## 3.1 頭份 PM2.5 2014-01-01 1 57
## 4.1 頭份 AMB_TEMP 2014-01-01 1 14
## 5.1 士林 NO2 2014-01-01 1 15
## 6.1 萬里 WIND_DIREC 2014-01-01 1 212
## 7.1 嘉義 PM10 2014-01-01 1 91
## 8.1 頭份 WIND_DIREC 2014-01-01 1 71
## 9.1 復興 WIND_SPEED 2014-01-01 1 .8
## 10.1 士林 WIND_SPEED 2014-01-01 1 1.2
## 11.1 二林 NOx 2014-01-01 1 23
## 12.1 士林 PM10 2014-01-01 1 56
## 13.1 復興 WIND_DIREC 2014-01-01 1 113
## 14.1 頭份 WIND_SPEED 2014-01-01 1 .9
## 15.1 士林 WIND_DIREC 2014-01-01 1 141
## 16.1 萬里 AMB_TEMP 2014-01-01 1 13
## 17.1 頭份 NO2 2014-01-01 1 24
## 18.1 復興 AMB_TEMP 2014-01-01 1 17
## 19.1 士林 AMB_TEMP 2014-01-01 1 16
## 20.1 頭份 NO 2014-01-01 1 2
## 21.1 士林 O3 2014-01-01 1 45
## 22.1 復興 PM2.5 2014-01-01 1 83
## 23.1 復興 RH 2014-01-01 1 74
## 24.1 士林 PM2.5 2014-01-01 1 28
## 25.1 前金 NO2 2014-01-01 1 30
## 26.1 頭份 NOx 2014-01-01 1 26
## 27.1 萬里 PM2.5 2014-01-01 1 30
## 28.1 士林 RH 2014-01-01 1 67
## 29.1 埔里 SO2 2014-01-01 1 2.5
## 30.1 二林 PM10 2014-01-01 1 78- 資料類型轉換
在寬版轉置前，需要進行2.2 資料察看的第2.點->「資料轉換」
首先，通常在資料整理時，必須參考資料來源單位所提供的說明檔。
在說明檔中有下列針對普通測站的說明：
#表示儀器檢核為無效值*表示程式檢核為無效值x表示人工檢核為無效值NR表示無降雨空白表示缺值風向資料888代表無風，999則代表儀器故障。
我可以利用grepl()來查看是否有特殊字串：
q1 = grepl("#", df.long$測量值, fixed=T)
sum(q1)## [1] 0q1 = grepl("*", df.long$測量值, fixed=T)
sum(q1)## [1] 0q1 = grepl("x", df.long$測量值, fixed=T)
sum(q1)## [1] 6720q1 = grepl("NR", df.long$測量值, fixed=T)
sum(q1)## [1] 0q1 = grepl("888", df.long$測量值, fixed=T)
sum(q1)## [1] 0q1 = grepl("999", df.long$測量值, fixed=T)
sum(q1)## [1] 0sum(nchar(df.long$測量值)==0)## [1] 0沒有異樣，可以進行轉換：
# 資料轉換
df.long[, "測量值old"] = df.long[, "測量值"]
df.long[, "測量值"] = suppressWarnings(as.numeric(df.long[, "測量值"]))
summary(df.long[, -c(1,2,3)])## 測量小時 測量值 測量值old
## Min. : 1.00 Min. : -1.00 Length:666144
## 1st Qu.: 6.75 1st Qu.: 3.20 Class :character
## Median :12.50 Median : 18.00 Mode :character
## Mean :12.50 Mean : 36.94
## 3rd Qu.:18.25 3rd Qu.: 50.00
## Max. :24.00 Max. :968.00
## NA's :9671# 日期資料轉換
df.long[, "測量日期"] = as.Date(df.long[, "測量日期"])
str(df.long)## 'data.frame': 666144 obs. of 6 variables:
## $ 測站名稱 : chr "士林" "士林" "頭份" "頭份" ...
## $ 測項英文名稱: chr "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## $ 測量日期 : Date, format: "2014-01-01" "2014-01-01" ...
## $ 測量小時 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ 測量值 : num 16 1.2 57 14 15 212 91 71 0.8 1.2 ...
## $ 測量值old : chr "16" "1.2" "57" "14" ...- 寬表(widw)：測站名稱，測量日期，測量小時，每個測量項目一欄
# 轉製成寬版
dfTemp = df.long
dfTemp[, "測量值old"] = NULL
idvar = c("測站名稱", "測量日期", "測量小時")
df.wide = reshape(dfTemp, idvar=idvar,
timevar="測項英文名稱",
v.names="測量值",
direction="wide")
str(df.wide)## 'data.frame': 56448 obs. of 15 variables:
## $ 測站名稱 : chr "士林" "頭份" "萬里" "嘉義" ...
## $ 測量日期 : Date, format: "2014-01-01" "2014-01-01" ...
## $ 測量小時 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ 測量值.NOx : num 16 26 5 37 75 23 32 17 46 20 ...
## $ 測量值.NO : num 1.2 2 0.9 3.1 29 1.4 1.3 1 5.8 2.5 ...
## $ 測量值.PM2.5 : num 28 57 30 65 83 46 85 69 59 24 ...
## $ 測量值.AMB_TEMP : num 16 14 13 15 17 13 16 14 17 15 ...
## $ 測量值.NO2 : num 15 24 4.1 34 46 21 30 16 41 17 ...
## $ 測量值.WIND_DIREC: num 141 71 212 61 113 38 34 312 128 118 ...
## $ 測量值.PM10 : num 56 76 38 91 134 78 120 79 76 40 ...
## $ 測量值.WIND_SPEED: num 1.2 0.9 4.7 1.2 0.8 1.8 1.7 0.9 1.5 0.9 ...
## $ 測量值.O3 : num 45 7.9 36 7.7 1.5 20 14 17 13 22 ...
## $ 測量值.RH : num 67 80 86 80 74 81 78 80 70 76 ...
## $ 測量值.SO2 : num 6.3 2.8 1.7 5.1 7.6 4.3 4.9 2.5 3.1 2 ...
## $ 測量值.CO : num 0.44 0.53 0.25 0.88 1.48 0.44 0.77 0.65 0.87 0.53 ...
## - attr(*, "reshapeWide")=List of 5
## ..$ v.names: chr "測量值"
## ..$ timevar: chr "測項英文名稱"
## ..$ idvar : chr [1:3] "測站名稱" "測量日期" "測量小時"
## ..$ times : chr [1:12] "NOx" "NO" "PM2.5" "AMB_TEMP" ...
## ..$ varying: chr [1, 1:12] "測量值.NOx" "測量值.NO" "測量值.PM2.5" "測量值.AMB_TEMP" ...# 整理欄位名稱
colnames(df.wide)[-c(1,2,3)] = substr(colnames(df.wide)[-c(1,2,3)],
5, nchar(colnames(df.wide)[-c(1,2,3)]))
attr(df.wide, "reshapeWide") = NULL
str(df.wide)## 'data.frame': 56448 obs. of 15 variables:
## $ 測站名稱 : chr "士林" "頭份" "萬里" "嘉義" ...
## $ 測量日期 : Date, format: "2014-01-01" "2014-01-01" ...
## $ 測量小時 : int 1 1 1 1 1 1 1 1 1 1 ...
## $ NOx : num 16 26 5 37 75 23 32 17 46 20 ...
## $ NO : num 1.2 2 0.9 3.1 29 1.4 1.3 1 5.8 2.5 ...
## $ PM2.5 : num 28 57 30 65 83 46 85 69 59 24 ...
## $ AMB_TEMP : num 16 14 13 15 17 13 16 14 17 15 ...
## $ NO2 : num 15 24 4.1 34 46 21 30 16 41 17 ...
## $ WIND_DIREC: num 141 71 212 61 113 38 34 312 128 118 ...
## $ PM10 : num 56 76 38 91 134 78 120 79 76 40 ...
## $ WIND_SPEED: num 1.2 0.9 4.7 1.2 0.8 1.8 1.7 0.9 1.5 0.9 ...
## $ O3 : num 45 7.9 36 7.7 1.5 20 14 17 13 22 ...
## $ RH : num 67 80 86 80 74 81 78 80 70 76 ...
## $ SO2 : num 6.3 2.8 1.7 5.1 7.6 4.3 4.9 2.5 3.1 2 ...
## $ CO : num 0.44 0.53 0.25 0.88 1.48 0.44 0.77 0.65 0.87 0.53 ...head(df.wide)## 測站名稱 測量日期 測量小時 NOx NO PM2.5 AMB_TEMP NO2 WIND_DIREC PM10
## 1.1 士林 2014-01-01 1 16 1.2 28 16 15.0 141 56
## 3.1 頭份 2014-01-01 1 26 2.0 57 14 24.0 71 76
## 6.1 萬里 2014-01-01 1 5 0.9 30 13 4.1 212 38
## 7.1 嘉義 2014-01-01 1 37 3.1 65 15 34.0 61 91
## 9.1 復興 2014-01-01 1 75 29.0 83 17 46.0 113 134
## 11.1 二林 2014-01-01 1 23 1.4 46 13 21.0 38 78
## WIND_SPEED O3 RH SO2 CO
## 1.1 1.2 45.0 67 6.3 0.44
## 3.1 0.9 7.9 80 2.8 0.53
## 6.1 4.7 36.0 86 1.7 0.25
## 7.1 1.2 7.7 80 5.1 0.88
## 9.1 0.8 1.5 74 7.6 1.48
## 11.1 1.8 20.0 81 4.3 0.442.7 儲存檔案
# 輸出檔案名稱
outDir = "./data/abc"
p1 = gregexpr(".csv", filename)
# 利用檔名製作輸出檔名
outName = paste(outDir, "/L",
substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
".rds", sep="")
saveRDS(df.long, outName)
# 輸出檔案名稱
outDir = "./data/abc"
p1 = gregexpr(".csv", filename)
# 利用檔名製作輸出檔名
outName = paste(outDir, "/W",
substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
".rds", sep="")
saveRDS(df.wide, outName)這樣檔案就輸出至abc資料夾裡了：
2.8 寫成函數(並大量作業)
由於有眾多檔案，我們不可能一個一個整理輸出，此時我們將上述程式碼整理成函數，使該函數可以大量作業，無須我們擠一個一個慢慢來。
程式碼如下：
pre1 = function(filename, outDir){
# 載入套件
library(dplyr)
# 讀入資料
df = read.csv(filename, fileEncoding="UTF-8-BOM")
if (sum(nchar(df$X.monitordate.)!=10)) df$X.monitordate. = substr(df$X.monitordate., 1, 10)
# 防止重複資料
df = unique(df)
# 測站資訊檔
stateInfo = data.frame(df$X.siteid., df$X.sitename.)
colnames(stateInfo) = c("測站代碼", "測站名稱")
stateInfo = unique(stateInfo)
# 測項資訊檔
itemInfo = data.frame(df$X.itemid., df$X.itemname.,
df$X.itemengname., df$X.itemunit.)
colnames(itemInfo) = c("測項代碼", "測項中文名稱", "測項英文名稱", "測項單位")
itemInfo = unique(itemInfo)
# 更改欄位名稱
dataInfo = dplyr::select(
df, X.sitename., X.itemengname.,
X.monitordate., starts_with("X.monitorvalue")
)
colnames(dataInfo) = c("測站名稱", "測項英文名稱", "測量日期",
paste("h0",seq(0,9),sep=""),
paste("h",seq(10,23),sep=""))
# 轉製成長版
varying = c(paste("h0",seq(0,9),sep=""),
paste("h",seq(10,23),sep=""))
df.long = reshape(dataInfo, varying=varying,
v.names="測量值",
timevar="測量小時",
direction="long")
df.long$id = NULL
attr(df.long, "reshapeLong") = NULL
df.long[, "測量值old"] = df.long[, "測量值"]
df.long[, "測量值"] = suppressWarnings(as.numeric(df.long[, "測量值"]))
# 日期資料轉換
df.long[, "測量日期"] = as.Date(df.long[, "測量日期"])
# 排序
# 轉製成寬版
dfTemp = df.long
dfTemp[, "測量值old"] = NULL
idvar = c("測站名稱", "測量日期", "測量小時")
df.wide = reshape(dfTemp, idvar=idvar,
timevar="測項英文名稱",
v.names="測量值",
direction="wide")
# 整理欄位名稱
colnames(df.wide)[-c(1,2,3)] = substr(colnames(df.wide)[-c(1,2,3)],
5, nchar(colnames(df.wide)[-c(1,2,3)]))
attr(df.wide, "reshapeWide") = NULL
#===長板===
# 利用檔名製作輸出檔名
p1 = gregexpr(".csv", filename)
outName = paste(outDir, "/L",
substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
".rds", sep="")
saveRDS(df.long, outName)
#===寬板===
# 輸出檔案名稱
p1 = gregexpr(".csv", filename)
# 利用檔名製作輸出檔名
outName = paste(outDir, "/W",
substr(filename, p1[[1]][1]-8, p1[[1]][1]-2),
".rds", sep="")
saveRDS(df.wide, outName)
# 輸出
return(list(stateInfo=stateInfo, itemInfo=itemInfo))
}測試 pre1() 函數：
filename = "./data/G03_P13_2014_2022/空氣品質監測小時值資料(一般污染物,每日更新) (2014-01).csv"
outDir="./data/S_0165_P13_G03"
q1 = pre1(filename, outDir=outDir)
head(q1$stateInfo)## 測站代碼 測站名稱
## 1 11 士林
## 3 25 頭份
## 6 3 萬里
## 7 42 嘉義
## 9 71 復興
## 11 35 二林q1$itemInfo## 測項代碼 測項中文名稱 測項英文名稱 測項單位
## 1 5 氮氧化物 NOx ppb
## 2 6 一氧化氮 NO ppb
## 3 33 細懸浮微粒 PM2.5 μg/m3
## 4 14 溫度 AMB_TEMP ℃
## 5 7 二氧化氮 NO2 ppb
## 6 11 風向 WIND_DIREC degrees
## 7 4 懸浮微粒 PM10 μg/m3
## 9 10 風速 WIND_SPEED m/sec
## 21 3 臭氧 O3 ppb
## 23 38 相對濕度 RH percent
## 29 1 二氧化硫 SO2 ppb
## 31 2 一氧化碳 CO ppm利用
list.dirs()查詢出某一目錄下的全部目錄名稱。- 查出來的第一個目錄是本身。
利用
list.files()查詢出某一目錄下的全部檔案名稱。
接著，就可以進行2014-2022的大量轉置作業~
補充：由於共9年，每年12個月，所以要輸入的數字為[1:108] !
# 檔案所在目錄
target = "./data/G03_P13_2014_2022/"
stateTemp = NULL
itemTemp = NULL
outDir = "./data/S_0165_P13_G03"
q2 = list.files(target, full.names=TRUE)
for (fName in q2[8]){
temp = tryCatch(
pre1(fName, outDir=outDir),
warning=function(w) print(paste(w,fName))
)
stateTemp = unique(rbind(stateTemp, temp$stateInfo))
itemTemp = unique(rbind(itemTemp, temp$itemInfo))
}
saveRDS(stateTemp, "./data/state_p13.rds")
2.9 年度資料
上一章節(2.8)的資料是每一年度每一個月的所有測站資料。
此章節主要是要將資料結合成每一年度所有測站的資料。
針對長版資料：
因為欄位固定，各檔案的欄位均相同。
將各檔案讀入後，再利用
rbind()指令合併即可。
針對寬版資料：
因各測站的測項並不固定，所以寬版檔案欄位不固定。
將各檔案讀入後，再利用
dplyr::bind_rows()指令合併。
先建立資料夾儲存檔案：
S_0165_p_13_year
程式碼如下：
setwd("E:/rumiworkspace")
# 設定輸入資料夾及輸出資料夾
target = "./data/S_0165_P13_G03/"
outDir = "./data/S_0165_p_13_year/"
# 設定長版或寬版
type="W"
fileList = list.files(target, full.names=TRUE)
q1 = unlist(gregexpr(".rds", fileList))
t1 = substr(fileList,q1-8,q1-8)
t2 = substr(fileList,q1-7,q1-4)
yearList = unique(t2[t1==type])
for (year in unique(yearList)){
fileListTemp = fileList[(t1==type) & (t2==year)]
temp = NULL
for (filename in fileListTemp){
# 讀入檔案
# 長版檔案使用
# 寬版檔案使用
temp = dplyr::bind_rows(temp, readRDS(filename))
}
# 寫出
# 組合檔名
outFile = paste(outDir, type, year, ".rds", sep="")
saveRDS(temp, outFile)
}
2.10 補充事項
在進行以上9個章節內容時，可能會發生的錯誤和解決辦法!
(這邊只提較多人出現的情況)
- Error in pre1(filename, outDir = outDir) : could not find function “pre1”
Ans：
未創建此函數、未執行此函數、函數名命名與執行寫入函數名不一致。
Error in file(file, “rt”, encoding = fileEncoding) :
cannot open the connection
Ans：
大多人都是在filename、outDir路徑上出現問題，會去檢查你的資料路徑和名稱是否與你寫入程式一致。
- 錯誤: The name of the input file cannot contain the special shell characters: [ <>()|\:&;#?*’] (attempted to copy to a version without those characters ‘Oscar-s-analysis-report.Rmd’ however that file already exists) 停止執行
Ans：
rmarkdown將刪除這些字符以嘗試沒有這些字符的另一個文件。如果我沒記錯的話，我們需要它來正確處理資源。回到儲存該檔案資料夾，將錯誤且多餘的檔案並重新開啟即可。
- 以些人在跑大量函數時會遇到，明明其他資料都很順的跑完了，唯獨零星幾個資料跑不出來。
Ans：
我的處理辦法是檢查資料名稱、路徑與寫入程式的是否一致。
若一致，則我認為是資料異常，重新下載資料再執行一次即可。
3.資料庫建置
完成第二章節內容後，我將要把寫好的資料寫入資料庫裡。
3.1 載入RMySQL
library(RMySQL)## Warning: 套件 'RMySQL' 是用 R 版本 4.2.3 來建造的## 載入需要的套件：DBI## Warning: 套件 'DBI' 是用 R 版本 4.2.3 來建造的3.2 建立資料表
設定連接管道，並將此管道命名為`con`
設定編碼，以避免出現亂碼
設定工作資料庫
才能將資料放入在我們指定的資料庫。
製作表格名稱
注意：將 - 換成 _底線。
- 找出適當的SQL語法
| 可以利用phpmyadmin中的SQL視窗試執行。 |
|---|
CREATE TABLE IF NOT EXISTS 表格名稱( 測站名稱 VARCHAR(3), 測項英文名稱 VARCHAR(15), 測量日期 DATE, 測量小時 TINYINT, 測量值 FLOAT, 測量值old VARCHAR(10), PRIMARY KEY (測站名稱, 測量日期, 測量小時) ); |
再利用套件RMySQL的適當函數執行。
資料插入資料表
關閉連機管道(這個步驟很重要)
否則會創建過多管道，使之後要執行程式會出現下面這錯誤!
如何解決錯誤：
Error in .local(drv, …) : Cannot allocate a new connection: 16 connections already opened
解決辦法：
library(RMySQL)
cons <- dbListConnections(MySQL())
for(con in cons)dbDisconnect(con)建立資料表的程式碼如下：
userName = "asis"
pswd = "12345678"
hostID = "127.0.0.1"
#source("./myFun/sqlConnectMac.R", local=TRUE)
con = dbConnect(RMySQL::MySQL(),
user=userName,
password=pswd,
host=hostID)
# 設定編碼
strQuery = "SET NAMES UTF8;"
#strQuery = "SET NAMES big5;" # 學校電腦教室需使用 big5
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 設定工作資料庫
# 若無此資料庫，請先建立
strQuery = "USE airp;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 讀取資料
setwd("E:/rumiworkspace")
filename = "./data/S_0165_P13_G03/L2014-05.rds"
df = readRDS(filename)
# 製做表格名稱
t1 = unlist(gregexpr(".rds", filename))
tbname = substr(filename, t1-8, t1-1)
# 將 - 換成 _底線
# 資料庫系統中，無法接受『-』為檔案名稱
tbname = gsub("-", "_", tbname)
# 設定查詢語法
strQuery = sprintf("
CREATE TABLE IF NOT EXISTS %s (
測站名稱 VARCHAR(3),
測項英文名稱 VARCHAR(15),
測量日期 DATE,
測量小時 TINYINT,
測量值 FLOAT,
測量值old VARCHAR(10),
PRIMARY KEY (測站名稱, 測量日期, 測量小時)
);", tbname
)
strQuery## [1] "\nCREATE TABLE IF NOT EXISTS L2014_05 (\n 測站名稱 VARCHAR(3),\n 測項英文名稱 VARCHAR(15),\n 測量日期 DATE,\n 測量小時 TINYINT,\n 測量值 FLOAT,\n 測量值old VARCHAR(10),\n PRIMARY KEY (測站名稱, 測量日期, 測量小時)\n);"# 執行查詢並回收執行結果
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 瀏覽查詢結果
#queryResult
# 關閉連結管道
queryResult = dbDisconnect(conn=con)3.3 dbWriteTable()
利用RMySQL套件的函數dbWriteTable()插入資料。
函數dbWriteTable()的功用是將R環境中的
資料框(data frame)物件
外部檔案
寫入MySQL資料庫成為資料表，語法如下：
{dbWriteTable(conn, name, value,} field.types=NULL, row.names=TRUE, overwrite=FALSE, append=FALSE, header=TRUE, sep=",", ...)
參數：
conn：一個有效的MySQL數據庫連接對象，由dbConnect函數創建。name：要創建的新表的名稱。value：要寫入表中的資料框（data frame）。row.names：一個邏輯值，表示是否將資料框的行名作為MySQL表的一個列。預設值為TRUE，表示將行名作為一個列。field.types：一個可選的參數，用於指定要在MySQL表中創建的各列的數據類型。它是一個命名字符向量，其中鍵是資料框的列名，值是要創建的MySQL列的數據類型。
返回值：
TRUE表示成功將資料框寫入到MySQL表中；
FALSE表示寫入失敗。
3.4 寫成函數
接著我們來接續上面的動作，並寫成函數執行。
tbCreate = function(filename){
# 載入套件
library(RMySQL)
# 設定連接管道，並將此管道命名為`con`
# 注意：帳秘每個人都不一樣
# host 的設定也有可能不一樣
userName = "asis"
pswd = "12345678"
hostID = "127.0.0.1"
#source("./myFun/sqlConnectMac.R", local=TRUE)
con = dbConnect(RMySQL::MySQL(),
user=userName,
password=pswd,
host=hostID)
# 設定編碼
strQuery = "SET NAMES UTF8;"
#strQuery = "SET NAMES big5;" # 學校電腦教室需使用 big5
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 設定工作資料庫
# 若無此資料庫，請先建立
strQuery = "USE airp;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 讀取資料
#filename = "./data/aqx_p_13/L2014-01.rds"
df = readRDS(filename)
# 製做表格名稱
t1 = unlist(gregexpr(".rds", filename))
tbname = substr(filename, t1-8, t1-1)
# 將 - 換成 _底線
tbname = gsub("-", "_", tbname)
# 資料插入資料表
bWriteTable(
conn=con,
name=tbname,
value=df,
row.names=FALSE,
overwrite=TRUE,
append=FALSE
)
# 關閉連結管道
queryResult = dbDisconnect(conn=con)
}3.5 大量作業
程式碼與2.8 寫成函數(並大量作業)此一節相似。
註：這邊先以2014年1月作示範，其餘皆有另行下載。
#setwd("E:/rumiworkspace")
#檔案所在目錄(每個人路徑資料夾名稱不同!!!)
#target = "./data/S_0165_P13_G03"
#查出目錄內所有檔案名稱
#q2 = list.files(target, full.names=TRUE)
#找出長版檔案名稱
#type = "L"
#fileList = q2[grepl(type, q2)]
# for (fName in fileList[1]){
# temp = tryCatch(
# tbCreate(fName),
# warning=function(w) print(paste(w,fName))
# )
# }
#}可能發生的錯誤：
- Error in .local(conn, statement, …) : could not run statement: Table ‘l2014_06’ already exists
解決辦法：資料來源路徑出現問題，修改路徑
- cannot open file ‘./myFun/sqlConnectMac.R’: No such file or directoryError in file(filename, “r”, encoding = encoding) : cannot open the connection
解決辦法：#source(“./myFun/sqlConnectMac.R”, local=TRUE)
3.6 補充說明
R中與MySQL數據庫進行連接並操作的優點：
通過與MySQL連接，可以利用R的強大功能對數據進行清洗、轉換、分析和可視化。
MySQL是一個廣泛使用的關聯式數據庫管理系統，具有高效的查詢和數據存儲能力。通過將數據輸入到MySQL中，可以利用MySQL的高效能和靈活性進行數據管理、查詢和後續分析。
將數據存儲在MySQL數據庫中，可以實現數據的共享和協作。多個使用者可以通過連接到MySQL數據庫，共同訪問和分析數據，從而促進團隊間的合作。
例如：陳明輝老師分享給11102-17201統計資料庫班上ASIS資料庫，幫忙學生們都整理好資料，以供我們在[4.研究主題]上好使用。
是否除了RMySQL套件，還有其他一些套件可以實現在R中與MySQL數據庫進行連接並輸入資料的功能或類似之?
RMariaDB：這是另一個在R中與MySQL數據庫交互的套件，提供類似於RMySQL的功能。示範： library(RMariaDB)
# 建立與MySQL數據庫的連接
con <- dbConnect(RMariaDB::MariaDB(), user = “username”, password = “password”,
dbname = “database_name”, host = “localhost”)
# 創建一個資料框
data <- data.frame(ID = c(1, 2, 3), Name = c(“Alice”, “Bob”, “Charlie”))
# 將資料框寫入到MySQL數據庫中的新表
dbWriteTable(con, name = “my_table”, value = data)
# 斷開與MySQL數據庫的連接
dbDisconnect(con)
RSQLite：這是一個用於在R中與SQLite數據庫進行交互的套件，它提供了類似於RMySQL的函數，如dbWriteTable可以將資料框寫入到SQLite數據庫。DBI套件：這是一個通用的數據庫接口套件，它提供了一致的API用於在R中與各種數據庫進行交互，包括MySQL、SQLite、PostgreSQL等。通過使用DBI套件，可以使用通用的函數來操作不同的數據庫。
有時在執行 ” library(RMySQL) “ 程式時會出現此提醒
載入需要的套件：DBI
這些套件都提供了類似的函數和功能，可以根據具體的需求和數據庫系統選擇適合的套件。
4.研究主題
4.1 摘要
本研究以環保署空氣品質監測之資料進行試算，討論空氣品質之變化情形。
何謂空氣品質?
空氣品質是指大氣中的空氣成分和污染物的特性和濃度，以及對人類健康和環境的影響程度。
國際上常用的空氣品質指標包括空氣質量指數（Air Quality Index，AQI），就是整合過去空氣污染指標(PSI) 係由每日監測所得懸浮微粒、二氧化硫、一氧化碳、臭氧及二氧化氮等5種主要污染物之濃度值換算而得之空氣污染指標值。 及細懸浮微粒(PM2.5)指標。
AQI值0-50代表空氣品質良好(綠色)、51-100為普通(黃色)、101-150為對敏感族群不良(橘色)、151-200為對所有族群不良(紅色)、201-300為非常不良(紫色)、301-500為有害(褐紅)。從而提供一個簡單明瞭的評估標準，幫助人們了解和應對空氣污染問題。
4.2 研究動機與目的
藉由數據資料，深入了解台灣的空氣品質現狀(所以此主題抓取2022的資料作研究)、污染源以及對人體健康和環境的影響。再以不同方面去探討(所以分成兩大研究主題)，通過分析不同地區和城市的空氣質量數據，我們可以評估區域情況並給予建議。
4.3 研究主題一
探討 汙染物的季節變化趨勢 及 地區空氣品質的影響因素分析
4.3.1 子主題一
不同季節下PM2.5和PM10濃度的變化趨勢及分析建議
4.3.2 研究主題說明
根據葉光芃在2019年的10月發表的最新論文指出 ，“粗顆粒帶來的健康效應鮮少受到關注，PM 2.5-10雖無法直達人體內肺泡，卻可滯留在上呼吸道內，會導致發炎和細胞中毒反應，可預見將導致越來越多死亡率，以及呼吸道與心血管疾病的就診率”。
所以此子主題針對pm2.5和pm10這兩測項資料，來看不同季節的變化趨勢，再以不同側站(地區)哪裡是最嚴重(濃度最高)的，深入討論該地區特徵，來得出結論。
4.3.3 方法及結果說明
連接數據庫
userName = "s10170165"
pswd = "s10170165@MCU1234"
hostID = "120.125.72.30"
con = dbConnect(RMySQL::MySQL(),
user=userName,
password=pswd,
host=hostID)
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
#選擇要用數據庫
dbSendQuery(con, "USE airp15")## <MySQLResult:878929267,1,1>獲取日的平均PM2.5和PM10濃度數據(柱狀圖)
SELECT 測量日期, AVG(`PM2.5`), AVG(`PM10`)
FROM W2022_dailyMean
GROUP BY 測量日期
;| 測量日期 | AVG(PM2.5) |
AVG(PM10) |
|---|---|---|
| 2022-01-01 | 15.46741 | 27.14263 |
| 2022-01-02 | 16.43269 | 29.28152 |
| 2022-01-03 | 25.44789 | 43.46165 |
| 2022-01-04 | 20.65050 | 34.35375 |
| 2022-01-05 | 24.08404 | 38.46897 |
| 2022-01-06 | 16.69693 | 30.79240 |
| 2022-01-07 | 13.95803 | 26.29919 |
| 2022-01-08 | 15.28793 | 29.62981 |
| 2022-01-09 | 16.50811 | 30.11083 |
| 2022-01-10 | 15.38953 | 27.45594 |
儲存結果在R中數據框
# 執行SQL查詢
datap1 <- dbGetQuery(con,
"SELECT 測量日期, AVG(`PM2.5`) 'PM2.5' , AVG(`PM10`) 'PM10'
FROM W2022_dailyMean
GROUP BY 測量日期;")
# 打印结果
head(datap1)## 測量日期 PM2.5 PM10
## 1 2022-01-01 15.46741 27.14263
## 2 2022-01-02 16.43269 29.28152
## 3 2022-01-03 25.44789 43.46165
## 4 2022-01-04 20.65050 34.35375
## 5 2022-01-05 24.08404 38.46897
## 6 2022-01-06 16.69693 30.79240柱狀圖繪圖
全台日平均PM2.5柱狀圖
library(ggplot2)## Warning: 套件 'ggplot2' 是用 R 版本 4.2.3 來建造的# 繪製PM2.5柱狀圖
ggplot(datap1, aes(x = 測量日期, y = PM2.5)) +
geom_bar(stat = "identity", fill = "blue", width = 0.5) +
labs(x = "測量日期", y = "平均PM2.5濃度") +
theme_minimal()
全台日平均PM10柱狀圖
library(ggplot2)
# 繪製PM10柱狀圖
ggplot(datap1, aes(x = 測量日期, y = PM10)) +
geom_bar(stat = "identity", fill = "green", width = 0.5) +
labs(x = "測量日期", y = "平均PM10濃度") +
theme_minimal()
柱狀圖分析
依分析後兩張柱狀圖配合這張全臺細懸浮微粒 (PM2.5)
標準監測數據平均四季變化圖(來源：空氣品質監測網)及PM10歷年月份平均值監測數據變化圖。
可以發現夏秋兩季全臺灣污染物平均濃度較低，而冬季及春季平均濃度較高。
獲取季節級別的數據(折線圖)
數據庫中只有日、小時級別的數據而沒有季節級別的數據，我們可以使用SQL查詢和計算來從日級別數據中推導季節級別的數據。以下展示如何計算季節級別的數據：
SELECT
YEAR(`測量日期`) AS 年份,
QUARTER(`測量日期`) AS 季度,
AVG(`PM2.5`) AS 平均PM2_5濃度,
AVG(`PM10`) AS 平均PM10濃度
FROM
W2022_dailyMean
GROUP BY
YEAR(測量日期),
QUARTER(測量日期)
ORDER BY
年份,
季度;| 年份 | 季度 | 平均PM2_5濃度 | 平均PM10濃度 |
|---|---|---|---|
| 2022 | 1 | 18.193016 | 31.16467 |
| 2022 | 2 | 11.274614 | 22.96019 |
| 2022 | 3 | 9.930115 | 22.55441 |
| 2022 | 4 | 13.119943 | 29.31640 |
上述代碼使用了YEAR()和QUARTER()函數從日期中提取年份和季度信息，併計算了每個季度的平均PM2.5和PM10濃度。
儲存結果在R中數據框
datap2 <- dbGetQuery(con,
"SELECT
YEAR(`測量日期`) AS 年份,
QUARTER(`測量日期`) AS 季度,
AVG(`PM2.5`) AS 平均PM2_5濃度,
AVG(`PM10`) AS 平均PM10濃度
FROM
W2022_dailyMean
GROUP BY
YEAR(測量日期),
QUARTER(測量日期)
ORDER BY
年份,
季度;")折線圖繪圖
全台季節PM2.5變化趨勢圖
library(ggplot2)
# 讀取從SQL查詢得到的季節級別數據
seasonal_data <- datap2
# 繪製季節級別的PM2.5變化趨勢圖
ggplot(seasonal_data, aes(x = 季度, y = 平均PM2_5濃度)) +
geom_line(color = "blue") +
geom_point(color = "blue") +
labs(x = "2022季度", y = "平均PM2.5濃度") +
theme_minimal()
全台季節PM10變化趨勢圖
library(ggplot2)
# 讀取從SQL查詢得到的季節級別數據
seasonal_data <- datap2
# 繪製季節級別的PM10變化趨勢圖
ggplot(seasonal_data, aes(x = 季度, y = 平均PM10濃度)) +
geom_line(color = "green") +
geom_point(color = "green") +
labs(x = "2022季度", y = "平均PM10濃度") +
theme_minimal()
全台季節PM2.5&PM10變化趨勢圖
ggplot(seasonal_data, aes(x = 季度)) +
geom_line(aes(y = 平均PM2_5濃度, color = "2.5")) +
geom_line(aes(y = 平均PM10濃度, color = "PM10")) +
labs(x = "Date", y = "Value") +
theme_minimal()
折線圖分析
此折線圖也可以發現夏秋兩季全臺灣污染物平均濃度較低，而冬季及春季平均濃度較高。
而PM10濃度比PM2.5濃度高。
獲取不同站點中PM2.5和PM10濃度最高值(熱力圖)
SELECT
`測站名稱`,
MAX(`PM2.5`) AS `最高PM2.5`,
MAX(`PM10`) AS `最高PM10`
FROM
`W2022_dailyMean`
GROUP BY
`測站名稱`
LIMIT 0, 25;| 測站名稱 | 最高PM2.5 | 最高PM10 |
|---|---|---|
| 三義 | 33.2917 | 51.5652 |
| 三重 | 39.2174 | 57.1739 |
| 中壢 | 43.2609 | 63.6957 |
| 中山 | 36.1667 | 52.7917 |
| 二林 | 56.9167 | 131.4350 |
| 仁武 | 58.9091 | 80.9167 |
| 冬山 | 32.5000 | 49.1667 |
| 前金 | 51.9583 | 81.6667 |
| 前鎮 | 59.1111 | 95.7222 |
| 南投 | 53.4167 | 84.2174 |
儲存結果在R中數據框
datap3 <- dbGetQuery(con,
"SELECT
`測站名稱`,
MAX(`PM2.5`) AS `最高PM2.5`,
MAX(`PM10`) AS `最高PM10`
FROM
`W2022_dailyMean`
GROUP BY
`測站名稱`
LIMIT 0, 25;")繪圖
PM2.5區域圖
PM2.5排名前五的地區：
嘉義、大寮、前鎮、仁武和善化
# 按照最高PM2.5排序，選擇排名前五的地區
top5 <- head(datap3[order(datap3$最高PM2.5, decreasing = TRUE), ], 5)
# 繪製板塊區域圖
ggplot(top5, aes(x = reorder(測站名稱, -最高PM2.5), y = 1, fill = 最高PM2.5)) +
geom_tile() +
labs(x = "測站名稱", y = NULL, fill = "最高PM2.5") +
scale_fill_gradient(low = "blue", high = "red") +
theme_minimal() +
theme(legend.position = "right")
PM10區域圖
PM10排名前五的地區： 二林、大寮、嘉義、前鎮和大里
# 按照最高PM10排序，選擇排名前五的地區
top5 <- head(datap3[order(datap3$最高PM10, decreasing = TRUE), ], 5)
# 繪製板塊區域圖
ggplot(top5, aes(x = reorder(測站名稱, -最高PM10), y = 1, fill = 最高PM10)) +
geom_tile() +
labs(x = "測站名稱", y = NULL, fill = "最高PM10") +
scale_fill_gradient(low = "blue", high = "red") +
theme_minimal() +
theme(legend.position = "right")
熱力圖程式碼解釋教學
使用 ggplot() 函數創建繪圖對象，並使用 geom_tile() 函數繪製熱力圖。在 aes() 函數中，將 測站名稱 作為x軸變量，使用1作為y軸變量（此處使用1是為了在圖表中創建一個單獨的行），並使用 最高PM2.5 作為填充顏色。通過 labs() 函數設置坐標軸和圖例的標籤，通過 scale_fill_gradient() 函數設置填充顏色的漸變範圍，最後使用 theme() 函數調整圖表的樣式。
熱力圖分析
根據熱力圖的呈現，可發現大寮、前鎮和嘉義這個測站皆重複於兩張圖上，表示這三個測站的地區在PM2.5和PM10上有些微嚴重，值得去探討其特徵。
大寮和前鎮 是在 高雄的測站
地區解釋(結論)
高雄
大寮和前鎮可能受到工業活動、交通污染和人口密度等因素的影響。
在這些地區，工廠排放、車輛尾氣和大量人口可能導致PM2.5和PM10濃度的升高。
進一步的研究可以關注這些測站所處地區的工業結構、交通流量以及人口分佈情況，以更詳細地了解造成高污染水平的具體原因。
嘉義
嘉義地區顯示出較高的PM2.5和PM10濃度，儘管該地區不屬於高雄地區，但嚴重性也很高。
在乾季(10 月至 4 月)嘉義位於東北季風隻被風處，細懸浮微粒較不容易擴散。
關閉連結通道
queryResult = dbDisconnect(conn=con)4.3.4 子主題二
各測站空氣汙染物比較
4.3.5 研究主題說明
以上針對完pm2.5和pm10，接著查看各測站汙染物占比，讓數據視覺化，以來看出各地區需特別注意的汙染物預防和趨勢。
這邊的各地區，我們分為北部、中部、雲嘉南、高屏和花東 空品區，再抓取以這五區的地方監測資料的測站綜合值來比較。
4.3.6 方法及結果說明
挑選測站
依據此網站地區中間測站和現有資料測站來挑選
北部：中正、板橋、觀音
中部：大里
雲嘉南：嘉義、崙背、埔里
高屏：鳳山、屏東
花東：花蓮、臺東
連接數據庫
userName = "s10170165"
pswd = "s10170165@MCU1234"
hostID = "120.125.72.30"
con = dbConnect(RMySQL::MySQL(),
user=userName,
password=pswd,
host=hostID)
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
#選擇要用數據庫
dbSendQuery(con, "USE airp15")## <MySQLResult:0,2,1>獲取各地區的地方監測站數據
北部
#北部
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('中正','板橋','觀音');| 測站名稱 | 測量日期 | AMB_TEMP | CH4 | CO | NMHC | NO | NO2 | NOx | O3 | PM10 | PM2.5 | RAINFALL | RH | SO2 | THC | WD_HR | WIND_DIREC | WIND_SPEED | WS_HR | PH_RAIN | RAIN_COND | RAIN_INT | CO2 | Wx.HR | Wy.HR | Wx.DIRECT | Wy.DIRECT |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 板橋 | 2022-01-01 | 17.2375 | 2.04375 | 0.398333 | 0.0875000 | 1.71250 | 15.2125 | 16.9833 | 29.0333 | 22.3750 | 12.00000 | 0.0000000 | 79.4583 | 1.520830 | 2.13125 | 120.2500 | 111.7500 | 1.04583 | 0.850000 | 0.00000 | 0.00000 | 0.0000000 | NA | 0.1142640 | 0.574366 | 0.0473125 | 0.745275 |
| 板橋 | 2022-01-02 | 17.3833 | 2.10542 | 0.430000 | 0.0979167 | 2.17917 | 16.2250 | 18.4500 | 24.4667 | 23.2917 | 12.45830 | 0.0000000 | 76.9167 | 0.812500 | 2.20333 | 96.0417 | 108.1670 | 1.27917 | 1.020830 | 0.00000 | 0.00000 | 0.0000000 | NA | -0.2328640 | 0.945135 | -0.2918940 | 1.132690 |
| 板橋 | 2022-01-03 | 18.8458 | 2.10000 | 0.483333 | 0.0758333 | 2.52500 | 16.3583 | 18.9083 | 38.6238 | 34.6250 | 20.00000 | 0.0000000 | 73.2083 | 0.545833 | 2.17583 | 97.8750 | 102.8750 | 1.53333 | 1.141670 | 0.00000 | 0.00000 | 0.0000000 | NA | -0.0833321 | 1.071400 | -0.1696980 | 1.349810 |
| 板橋 | 2022-01-04 | 18.5417 | 2.01792 | 0.340417 | 0.0425000 | 1.92500 | 13.5083 | 15.4792 | 36.4958 | 16.1250 | 9.29167 | 0.0000000 | 81.7083 | 1.837500 | 2.05792 | 93.9167 | 94.0417 | 2.38750 | 1.845830 | 0.00000 | 0.00000 | 0.0000000 | NA | -0.0024818 | 1.709700 | 0.0431246 | 2.178940 |
| 板橋 | 2022-01-05 | 19.3958 | 2.11292 | 0.622083 | 0.2375000 | 9.31250 | 23.8542 | 33.2083 | 18.9292 | 24.4167 | 14.58330 | 0.0000000 | 82.7083 | 1.470830 | 2.35042 | 195.5420 | 196.1250 | 1.16667 | 0.995833 | 0.00000 | 0.00000 | 0.0000000 | NA | 0.3928380 | -0.266902 | 0.4538090 | -0.209908 |
| 板橋 | 2022-01-06 | 17.1167 | 2.05316 | 0.415714 | 0.1036840 | 2.12105 | 17.6105 | 19.7842 | 26.9750 | 19.6522 | 10.34780 | 0.3636360 | 80.7917 | 0.933333 | 2.15684 | 101.0000 | 103.0870 | 2.15217 | 1.608700 | 1.29167 | 7.91250 | 0.3333330 | NA | -0.3617850 | 1.410720 | -0.4510510 | 1.846590 |
| 板橋 | 2022-01-07 | 16.2083 | 2.01625 | 0.378696 | 0.0941667 | 2.10000 | 18.2810 | 20.4143 | 28.0435 | 25.3333 | 11.54170 | 0.1875000 | 78.0000 | 1.065220 | 2.11042 | 109.9580 | 109.0420 | 1.77500 | 1.433330 | 1.09583 | 2.04167 | 0.2291670 | NA | -0.4550400 | 1.322350 | -0.4832120 | 1.637410 |
| 板橋 | 2022-01-08 | 16.5042 | 2.01542 | 0.316667 | 0.0420833 | 1.64167 | 12.4833 | 14.1708 | 32.1333 | 24.7500 | 10.25000 | 0.0000000 | 74.5833 | 1.175000 | 2.05750 | 103.8330 | 103.7920 | 1.79583 | 1.454170 | 0.00000 | 0.00000 | 0.0000000 | NA | -0.2655810 | 1.382720 | -0.2864310 | 1.675800 |
| 板橋 | 2022-01-09 | 17.8417 | 2.02417 | 0.335000 | 0.0462500 | 1.25417 | 12.3750 | 13.6792 | 33.2542 | 21.9167 | 11.20830 | 0.0000000 | 75.6250 | 1.225000 | 2.07042 | 104.3750 | 112.3750 | 1.51667 | 1.233330 | 0.00000 | 0.00000 | 0.0000000 | NA | -0.1944710 | 1.147880 | -0.3486550 | 1.360960 |
| 板橋 | 2022-01-10 | 17.6833 | 2.01667 | 0.400000 | 0.0837500 | 2.56667 | 16.9542 | 19.5708 | 27.7000 | 16.9167 | 9.58333 | 0.0208333 | 79.7083 | 1.000000 | 2.10042 | 111.4170 | 117.5830 | 1.77083 | 1.408330 | 0.23750 | 2.63750 | 0.0208333 | NA | -0.0200464 | 1.124950 | -0.0166594 | 1.385970 |
中部
#中部
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('大里');| 測站名稱 | 測量日期 | AMB_TEMP | CH4 | CO | NMHC | NO | NO2 | NOx | O3 | PM10 | PM2.5 | RAINFALL | RH | SO2 | THC | WD_HR | WIND_DIREC | WIND_SPEED | WS_HR | PH_RAIN | RAIN_COND | RAIN_INT | CO2 | Wx.HR | Wy.HR | Wx.DIRECT | Wy.DIRECT |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 大里 | 2022-01-01 | 18.7458 | 2.23042 | 0.384583 | 0.1320830 | 2.42917 | 15.5500 | 18.0458 | 20.3625 | 29.0000 | 18.8333 | 0.0000000 | 65.7917 | 1.62917 | 2.36250 | 132.250 | 178.667 | 1.77500 | 1.404170 | NA | NA | NA | 444.375 | 1.114320 | 0.0006554 | 1.2305000 | -0.0070190 |
| 大里 | 2022-01-02 | 18.6250 | 2.28167 | 0.411250 | 0.1612500 | 3.62500 | 13.7875 | 17.4625 | 20.0583 | 25.8750 | 17.2083 | 0.0000000 | 69.3333 | 1.61250 | 2.44292 | 240.333 | 253.667 | 1.60417 | 1.237500 | NA | NA | NA | 450.125 | 0.658168 | -0.1892270 | 0.8695330 | -0.3629620 |
| 大里 | 2022-01-03 | 19.1167 | 2.27833 | 0.559167 | 0.2208330 | 8.81250 | 20.0000 | 28.8667 | 23.2000 | 43.4167 | 27.2917 | 0.0000000 | 70.3333 | 1.86667 | 2.49917 | 265.250 | 262.000 | 1.55417 | 1.225000 | NA | NA | NA | 458.417 | 0.695231 | -0.3152430 | 0.9234620 | -0.4058590 |
| 大里 | 2022-01-04 | 19.4833 | 2.35083 | 0.684583 | 0.4237500 | 13.53330 | 26.4750 | 40.0708 | 14.8125 | 51.5000 | 30.2083 | 0.0000000 | 74.9167 | 2.03750 | 2.77458 | 230.458 | 228.542 | 1.12917 | 0.770833 | NA | NA | NA | 468.583 | 0.157299 | -0.4042500 | 0.0710028 | -0.4770700 |
| 大里 | 2022-01-05 | 19.9083 | 2.34565 | 0.601739 | 0.3600000 | 12.93480 | 23.6696 | 36.6435 | 14.2000 | 43.0455 | 26.3913 | 0.0000000 | 77.3750 | 1.76522 | 2.70565 | 209.500 | 206.458 | 1.33333 | 0.916667 | NA | NA | NA | 464.043 | 0.511665 | -0.2164610 | 0.6339280 | -0.2513420 |
| 大里 | 2022-01-06 | 18.7208 | 2.30053 | 0.416250 | 0.2357890 | 4.81250 | 15.2708 | 20.1292 | 19.2583 | 25.7083 | 15.7083 | 0.1666670 | 76.0833 | 1.07083 | 2.53632 | 223.667 | 172.750 | 2.21250 | 1.700000 | NA | NA | NA | 447.917 | 1.352810 | -0.1111270 | 1.6725000 | -0.0492783 |
| 大里 | 2022-01-07 | 17.5625 | 2.09429 | 0.304167 | 0.1035710 | 3.37500 | 12.3208 | 15.7208 | 24.1250 | 14.2083 | 7.1250 | 0.1250000 | 76.7917 | 1.37500 | 2.19786 | 218.125 | 186.917 | 2.25833 | 1.737500 | NA | NA | NA | 435.583 | 1.664270 | -0.1720170 | 1.9106000 | -0.2389280 |
| 大里 | 2022-01-08 | 18.2750 | 2.13417 | 0.322083 | 0.0770833 | 3.48750 | 13.0292 | 16.5667 | 23.2208 | 19.6250 | 10.7500 | 0.0416667 | 71.7917 | 1.65000 | 2.21125 | 259.250 | 291.958 | 1.92917 | 1.495830 | NA | NA | NA | 436.125 | 1.389560 | -0.3476450 | 1.7333900 | -0.3704840 |
| 大里 | 2022-01-09 | 18.5583 | 2.23208 | 0.358333 | 0.1120830 | 5.41250 | 13.9750 | 19.4167 | 24.6542 | 30.5833 | 17.6667 | 0.0000000 | 74.4583 | 1.44583 | 2.34417 | 230.333 | 220.958 | 1.77500 | 1.441670 | NA | NA | NA | 444.792 | 0.977188 | -0.3839520 | 1.1281800 | -0.5592080 |
| 大里 | 2022-01-10 | 18.5375 | 2.27542 | 0.467500 | 0.1841670 | 9.29583 | 17.0875 | 26.4208 | 19.6208 | 35.2083 | 20.4167 | 0.0000000 | 76.3750 | 1.48750 | 2.45958 | 231.792 | 222.042 | 1.60000 | 1.166670 | NA | NA | NA | 452.708 | 0.563592 | -0.3770280 | 0.7574100 | -0.5830180 |
雲嘉南
#雲嘉南
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('嘉義','崙背','埔里');| 測站名稱 | 測量日期 | AMB_TEMP | CH4 | CO | NMHC | NO | NO2 | NOx | O3 | PM10 | PM2.5 | RAINFALL | RH | SO2 | THC | WD_HR | WIND_DIREC | WIND_SPEED | WS_HR | PH_RAIN | RAIN_COND | RAIN_INT | CO2 | Wx.HR | Wy.HR | Wx.DIRECT | Wy.DIRECT |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 嘉義 | 2022-01-01 | 18.7333 | 2.04500 | 0.337083 | 0.121250 | 1.95417 | 13.3625 | 15.3542 | 25.2167 | 38.0000 | 21.8333 | 0.0000000 | 68.9583 | 2.38333 | 2.16625 | 45.7083 | 34.0417 | 2.63333 | 2.26250 | NA | NA | NA | NA | 1.94367 | 0.9097920 | 2.14516 | 1.1544500 |
| 嘉義 | 2022-01-02 | 18.8125 | 2.09083 | 0.361250 | 0.125000 | 2.37083 | 12.2083 | 14.6208 | 24.4500 | 43.3333 | 25.7083 | 0.0000000 | 74.2917 | 2.58750 | 2.21583 | 81.6667 | 83.0417 | 1.97917 | 1.72500 | NA | NA | NA | NA | 1.57099 | 0.3733600 | 1.73073 | 0.4685170 |
| 嘉義 | 2022-01-03 | 19.2458 | 2.13130 | 0.479565 | 0.175217 | 2.83913 | 18.3391 | 21.2217 | 26.7625 | 57.2273 | 34.1304 | 0.0000000 | 72.5417 | 3.75217 | 2.30652 | 75.2083 | 74.4167 | 2.20417 | 1.92083 | NA | NA | NA | NA | 1.78946 | 0.4826990 | 1.94530 | 0.5731490 |
| 嘉義 | 2022-01-04 | 19.7583 | 2.13833 | 0.497500 | 0.227917 | 3.86250 | 19.3375 | 23.2750 | 25.8095 | 63.6250 | 40.6250 | 0.0000000 | 76.8750 | 3.57917 | 2.36625 | 143.2920 | 155.5830 | 1.60417 | 1.43333 | NA | NA | NA | NA | 1.23180 | -0.0004508 | 1.30561 | 0.0397694 |
| 嘉義 | 2022-01-05 | 19.8708 | 2.17042 | 0.515417 | 0.252083 | 4.41250 | 21.7958 | 26.2500 | 15.8417 | 66.0417 | 44.2500 | 0.0000000 | 83.0833 | 3.23750 | 2.42250 | 93.5417 | 111.1670 | 1.81667 | 1.52917 | NA | NA | NA | NA | 1.21388 | 0.4234570 | 1.47956 | 0.3906820 |
| 嘉義 | 2022-01-06 | 19.1625 | 2.05833 | 0.378750 | 0.186667 | 3.00000 | 18.1042 | 21.1583 | 19.0375 | 40.7500 | 23.4167 | 0.0000000 | 79.8333 | 2.78333 | 2.24500 | 57.5833 | 60.0833 | 2.42917 | 2.17500 | NA | NA | NA | NA | 2.07063 | 0.3788190 | 2.27950 | 0.5202170 |
| 嘉義 | 2022-01-07 | 17.3583 | 2.01833 | 0.327083 | 0.147917 | 2.77083 | 16.6333 | 19.4583 | 18.8625 | 26.4167 | 16.9167 | 0.0416667 | 84.4583 | 2.30833 | 2.16625 | 49.5000 | 78.4167 | 2.36667 | 2.07083 | NA | NA | NA | NA | 1.86465 | 0.6074310 | 1.85713 | 0.5065030 |
| 嘉義 | 2022-01-08 | 18.3875 | 2.03917 | 0.350417 | 0.131667 | 2.44167 | 15.3042 | 17.7792 | 23.7875 | 37.3750 | 21.2083 | 0.0000000 | 75.9583 | 3.17500 | 2.17083 | 91.6667 | 162.6670 | 2.04583 | 1.68750 | NA | NA | NA | NA | 1.56178 | 0.3776430 | 1.87500 | 0.2891700 |
| 嘉義 | 2022-01-09 | 18.7458 | 2.07208 | 0.341667 | 0.112500 | 2.27083 | 12.7083 | 15.0042 | 28.2417 | 39.4583 | 24.5417 | 0.0000000 | 76.1667 | 4.16250 | 2.18458 | 102.2500 | 103.7080 | 2.17917 | 2.05000 | NA | NA | NA | NA | 1.93371 | 0.2691310 | 1.98095 | 0.3424650 |
| 嘉義 | 2022-01-10 | 18.7333 | 2.07625 | 0.372083 | 0.118750 | 3.10833 | 13.0000 | 16.1708 | 21.9524 | 41.3913 | 24.4783 | 0.0000000 | 79.7083 | 2.64167 | 2.19500 | 106.0830 | 127.7500 | 1.76250 | 1.57083 | NA | NA | NA | NA | 1.38046 | 0.1503190 | 1.45045 | 0.2574620 |
高屏
#高屏
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('鳳山','屏東');| 測站名稱 | 測量日期 | AMB_TEMP | CH4 | CO | NMHC | NO | NO2 | NOx | O3 | PM10 | PM2.5 | RAINFALL | RH | SO2 | THC | WD_HR | WIND_DIREC | WIND_SPEED | WS_HR | PH_RAIN | RAIN_COND | RAIN_INT | CO2 | Wx.HR | Wy.HR | Wx.DIRECT | Wy.DIRECT |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 屏東 | 2022-01-01 | 21.0250 | 2.15917 | 0.427917 | 0.1058330 | 1.212500 | 13.6708 | 14.9292 | 29.7667 | 41.6250 | 25.1250 | 0.000 | 59.5417 | 0.829167 | 2.26500 | 254.000 | 219.542 | 1.48750 | 1.22083 | NA | NA | NA | NA | 0.941790 | -0.3061450 | 1.0974000 | -0.2966800 |
| 屏東 | 2022-01-02 | 21.4958 | 2.15625 | 0.372083 | 0.0820833 | 0.829167 | 11.6375 | 12.5167 | 33.9708 | 32.8333 | 18.4583 | 0.000 | 62.5417 | 1.466670 | 2.23833 | 236.875 | 262.125 | 1.57083 | 1.27083 | NA | NA | NA | NA | 0.553664 | -0.3778020 | 0.5518180 | -0.5781180 |
| 屏東 | 2022-01-03 | 21.8333 | 2.21125 | 0.435000 | 0.1116670 | 1.704170 | 14.4875 | 16.2542 | 35.5417 | 45.0000 | 26.7500 | 0.000 | 63.3750 | 1.433330 | 2.32292 | 207.375 | 176.167 | 1.47917 | 1.11667 | NA | NA | NA | NA | 0.708287 | -0.0860183 | 0.8019740 | -0.0812447 |
| 屏東 | 2022-01-04 | 22.0583 | 2.21417 | 0.531667 | 0.1362500 | 2.112500 | 18.1542 | 20.3167 | 42.7458 | 61.2500 | 37.0833 | 0.000 | 66.1250 | 2.787500 | 2.35042 | 167.500 | 190.917 | 1.43333 | 1.16250 | NA | NA | NA | NA | 0.803308 | -0.0967507 | 0.6945240 | -0.1538960 |
| 屏東 | 2022-01-05 | 23.0042 | 2.17625 | 0.424167 | 0.1108330 | 0.908333 | 14.0375 | 14.9958 | 33.7750 | 56.8333 | 37.2917 | 0.000 | 68.5833 | 1.391670 | 2.28708 | 278.375 | 223.083 | 1.46250 | 1.20833 | NA | NA | NA | NA | 1.054390 | -0.3793600 | 1.1864700 | -0.3762350 |
| 屏東 | 2022-01-06 | 21.5708 | 2.15875 | 0.523750 | 0.1820830 | 3.645830 | 19.4292 | 23.1208 | 23.2583 | 64.6250 | 39.5000 | 0.000 | 72.7917 | 2.816670 | 2.34083 | 244.167 | 275.083 | 1.45000 | 1.15417 | NA | NA | NA | NA | 0.574779 | -0.6172840 | 0.7549260 | -0.7416020 |
| 屏東 | 2022-01-07 | 20.7750 | 2.09455 | 0.359545 | 0.1140910 | 2.166670 | 15.8095 | 18.0143 | 19.7792 | 36.0000 | 22.9048 | 0.025 | 74.6667 | 0.913636 | 2.20864 | 288.792 | 249.208 | 1.33750 | 1.11250 | NA | NA | NA | NA | 0.742434 | -0.6652090 | 0.8584420 | -0.6643760 |
| 屏東 | 2022-01-08 | 21.4667 | 2.13625 | 0.438750 | 0.1425000 | 2.916670 | 16.5750 | 19.5458 | 35.9714 | 53.8333 | 32.6667 | 0.000 | 70.8333 | 2.258330 | 2.27875 | 237.125 | 248.833 | 1.52083 | 1.12083 | NA | NA | NA | NA | 0.156468 | -0.6939870 | 0.0881084 | -0.9453870 |
| 屏東 | 2022-01-09 | 21.3000 | 2.13375 | 0.347083 | 0.0787500 | 1.404170 | 10.5000 | 11.9542 | 37.1458 | 47.2500 | 28.6250 | 0.000 | 69.2500 | 1.675000 | 2.21250 | 260.875 | 282.958 | 1.66250 | 1.30000 | NA | NA | NA | NA | 0.831068 | -0.7423670 | 0.8319360 | -1.0849300 |
| 屏東 | 2022-01-10 | 20.5208 | 2.12542 | 0.322917 | 0.0545833 | 2.212500 | 10.2000 | 12.4583 | 33.7750 | 42.5000 | 25.6250 | 0.000 | 72.6250 | 1.237500 | 2.18000 | 290.375 | 276.750 | 1.84167 | 1.43750 | NA | NA | NA | NA | 1.215780 | -0.6249600 | 1.4704800 | -0.7544330 |
花東
#花東
SELECT * FROM `W2022_dailyMean` WHERE `測站名稱` IN ('花蓮','臺東');| 測站名稱 | 測量日期 | AMB_TEMP | CH4 | CO | NMHC | NO | NO2 | NOx | O3 | PM10 | PM2.5 | RAINFALL | RH | SO2 | THC | WD_HR | WIND_DIREC | WIND_SPEED | WS_HR | PH_RAIN | RAIN_COND | RAIN_INT | CO2 | Wx.HR | Wy.HR | Wx.DIRECT | Wy.DIRECT |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 臺東 | 2022-01-01 | 19.8333 | NA | 0.222917 | NA | 1.437500 | 4.28750 | 5.75417 | 33.6542 | 11.7083 | 7.75000 | 0.0000000 | 77.3333 | 1.90833 | NA | 219.000 | 232.583 | 1.99167 | 1.63333 | 0.000000 | 0.000000 | 0.0000000 | NA | 1.1244600 | -0.1366420 | 1.352220 | -0.1598950 |
| 臺東 | 2022-01-02 | 20.5167 | NA | 0.202083 | NA | 0.658333 | 3.47500 | 4.19167 | 31.2500 | 13.2917 | 7.16667 | 0.0000000 | 74.0833 | 1.48750 | NA | 210.542 | 238.583 | 2.05417 | 1.62083 | 0.000000 | 0.000000 | 0.0000000 | NA | 1.2885100 | -0.2382630 | 1.605980 | -0.3151000 |
| 臺東 | 2022-01-03 | 20.6417 | NA | 0.261250 | NA | 1.352170 | 3.82174 | 5.21739 | 34.9292 | 17.2609 | 9.04348 | 0.0000000 | 75.2083 | 1.25833 | NA | 223.708 | 238.625 | 2.20417 | 1.86250 | 0.000000 | 0.000000 | 0.0000000 | NA | 1.4070100 | -0.3565540 | 1.677900 | -0.3411060 |
| 臺東 | 2022-01-04 | 20.8125 | NA | 0.263478 | NA | 1.117390 | 4.93478 | 6.08261 | 27.7333 | 12.2174 | 6.39130 | 0.0454545 | 82.0833 | 1.22174 | NA | 248.917 | 237.917 | 1.40000 | 1.16250 | 0.233333 | 0.358333 | 0.0208333 | NA | 0.9126610 | -0.3620370 | 1.070610 | -0.3234910 |
| 臺東 | 2022-01-05 | 22.3417 | NA | 0.266250 | NA | 1.141670 | 5.87083 | 7.07500 | 24.1167 | 13.1250 | 4.83333 | 0.0208333 | 84.8333 | 1.20000 | NA | 261.708 | 236.917 | 1.37500 | 0.96250 | 0.000000 | 0.000000 | 0.0000000 | NA | 0.0660739 | -0.3280410 | 0.169317 | -0.3324130 |
| 臺東 | 2022-01-06 | 20.7292 | NA | 0.282083 | NA | 1.537500 | 5.00417 | 6.56667 | 34.8042 | 12.2500 | 6.62500 | 0.0416667 | 81.0417 | 1.45417 | NA | 210.750 | 238.125 | 2.12083 | 1.76250 | 0.245833 | 1.070830 | 0.0208333 | NA | 1.4885700 | -0.1746700 | 1.775390 | -0.2681930 |
| 臺東 | 2022-01-07 | 19.9792 | NA | 0.250833 | NA | 1.216670 | 4.52917 | 5.78750 | 36.2833 | 10.7917 | 5.00000 | 0.0208333 | 75.6250 | 1.50417 | NA | 249.625 | 248.917 | 2.15833 | 1.77083 | 0.258333 | 0.587500 | 0.0208333 | NA | 1.4142800 | -0.3922460 | 1.706100 | -0.4985300 |
| 臺東 | 2022-01-08 | 19.6833 | NA | 0.228333 | NA | 0.912500 | 4.09167 | 5.04167 | 36.7833 | 10.2917 | 5.91667 | 0.0000000 | 76.6667 | 1.55833 | NA | 249.042 | 236.000 | 2.19583 | 1.81250 | 0.000000 | 0.000000 | 0.0000000 | NA | 1.3981000 | -0.4700110 | 1.650000 | -0.4798470 |
| 臺東 | 2022-01-09 | 20.4500 | NA | 0.214167 | NA | 0.850000 | 3.72083 | 4.61667 | 36.1792 | 14.6667 | 6.58333 | 0.0000000 | 75.9167 | 1.40000 | NA | 213.208 | 226.542 | 1.88333 | 1.57500 | 0.000000 | 0.000000 | 0.0000000 | NA | 1.1995200 | -0.2624550 | 1.460300 | -0.3206360 |
| 臺東 | 2022-01-10 | 21.0458 | NA | 0.232917 | NA | 1.695830 | 3.60833 | 5.34583 | 33.9708 | 15.7500 | 5.79167 | 0.0000000 | 78.7083 | 1.46250 | NA | 216.625 | 230.458 | 2.12083 | 1.73333 | 0.000000 | 0.000000 | 0.0000000 | NA | 1.2490200 | 0.0381417 | 1.561820 | 0.0688152 |
儲存結果在R中數據框
# 從數據庫中獲取各地區的地方監測站數據
region_data1 <- dbGetQuery(con, "SELECT SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('中正','板橋','觀音') GROUP BY `測站名稱`")
region_data1 <- colSums(region_data1)
region_data2 <- dbGetQuery(con, "SELECT SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('大里') GROUP BY `測站名稱`")
region_data2 <- colSums(region_data2)
region_data3 <- dbGetQuery(con, "SELECT SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('嘉義','崙背','埔里') GROUP BY `測站名稱`")
region_data3 <- colSums(region_data3)
region_data4 <- dbGetQuery(con, "SELECT SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('鳳山','屏東') GROUP BY `測站名稱`")
region_data4 <- colSums(region_data4)
region_data5 <- dbGetQuery(con, "SELECT SUM(`PM2.5`) AS `PM2.5`, SUM(`PM10`) AS `PM10`, SUM(`SO2`) AS `SO2`, SUM(`NO2`) AS `NO2`, SUM(`CO`) AS `CO`, SUM(`O3`) AS `O3` FROM W2022_dailyMean WHERE `測站名稱` IN ('花蓮','臺東') GROUP BY `測站名稱` ;")
region_data5 <- colSums(region_data5)畫圖
北部
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data1,xlab = "北部測站",ylab = "污染物濃度")
中部
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data2,xlab = "中部測站",ylab = "污染物濃度")
雲嘉南
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data3,xlab = "雲嘉南測站",ylab = "污染物濃度")
高屏
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data4,xlab = "高屏測站",ylab = "污染物濃度")
花東
# 可視化各地區的污染物總和
library(ggplot2)
barplot(region_data5,xlab = "花東測站",ylab = "污染物濃度")
結論分析
根據可視化結果，我們可以觀察到所有地區的O3和PM10污染物總和相對很高，而PM2.5次高。這可能意味著台灣面臨著較高的臭氧和顆粒物污染問題。
這裡特別提到台北和花東地區測站的PM10也是第二高，但是相對其他地區也明顯的低，那麼我將會在研究主題二有更深入的探討。
高O3濃度：地區可能受到較高的光照和溫暖的氣候影響，這促進了大氣中氮氧化物和揮發性有機化合物的光化反應，產生了更多的臭氧。這種情況通常在夏季和高溫季節更為明顯。
高PM10濃度：地區可能受到交通污染、工業排放和城市化程度高的影響。車輛排放、工業廢氣和揚塵等因素可能導致空氣中顆粒物（PM10）的積累和增加。
這些解釋是基於對污染物總和可視化結果的觀察和常見的空氣污染來源分析。
關閉連結通道
queryResult = dbDisconnect(conn=con)4.4 研究主題二
都會區(這邊針對台北市)臭氧(O3)每日最高值或惡化情況之分析
4.4.1 子主題一
臭氧濃度與其他空氣汙染物相關性
4.4.2 研究主題說明
首先，我們以環保署所設之空氣品質監測站中的台北市測站為對象，選取2022年1～12月空氣品質測定報告之資料，利用迴歸分析，加以分析及探討。
4.4.3 方法及結果說明
此節為分析架構 與 回歸模型策略
>>取得並整理資料
1. 查詢 W2022_dailyMean 資料表中，有那些是台北市測站(利用排除重複)
程式碼
strQuery = "SELECT DISTINCT 測站名稱 FROM W2022_dailyMean;"
結果
## 測站名稱
## 1 三義
## 2 三重
## 3 中壢
## 4 中山
## 5 二林
## 6 仁武
## 7 冬山
## 8 前金
## 9 前鎮
## 10 南投
## 11 古亭
## 12 善化
## 13 嘉義
## 14 土城
## 15 埔里
## 16 基隆
## 17 士林
## 18 大同
## 19 大園
## 20 大城
## 21 大寮
## 22 大里
## 23 安南
## 24 宜蘭
## 25 富貴角
## 26 小港
## 27 屏東
## 28 崙背
## 29 左營
## 30 平鎮
## 31 彰化
## 32 復興
## 33 忠明
## 34 恆春
## 35 斗六
## 36 新店
## 37 新港
## 38 新營
## 39 新竹
## 40 新莊
## 41 朴子
## 42 松山
## 43 板橋
## 44 林口
## 45 林園
## 46 桃園
## 47 楠梓
## 48 橋頭
## 49 永和
## 50 汐止
## 51 沙鹿
## 52 淡水
## 53 湖口
## 54 潮州
## 55 竹山
## 56 竹東
## 57 線西
## 58 美濃
## 59 臺南
## 60 臺東
## 61 臺西
## 62 花蓮
## 63 苗栗
## 64 菜寮
## 65 萬華
## 66 萬里
## 67 西屯
## 68 觀音
## 69 豐原
## 70 金門
## 71 關山
## 72 陽明
## 73 頭份
## 74 馬公
## 75 馬祖
## 76 鳳山
## 77 麥寮
## 78 龍潭再根據空氣品質監測網的測站資訊 和 我們排除重複得出的測站資訊，已得到我們之後所針對的測站有：
古亭、士林、大同、陽明、松山、中山、萬華
2. 抓取這些測站的測項資料，並製作成資料框。
程式碼
strQuery = "SELECT `O3`,`CH4`,`CO`,`NMHC`,`SO2`,`PM10`,`PM2.5`,`NO2`,`NO2`,`NO`,`THC`,`CO2`,`CH4` FROM W2022_dailyMean WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`CH4`IS NOT NULL
AND`CO`IS NOT NULL
AND`NMHC`IS NOT NULL
AND`SO2`IS NOT NULL
AND`PM10`IS NOT NULL
AND`PM2.5`IS NOT NULL
AND`NO2`IS NOT NULL
AND`NO2`IS NOT NULL
AND`NO`IS NOT NULL
AND`THC`IS NOT NULL
AND`CO2`IS NOT NULL
AND`CH4`IS NOT NULL;"
結果
## O3 CH4 CO NMHC SO2 PM10 PM2.5 NO2 NO
## 1 20.3625 2.23042 0.384583 0.132083 1.62917 29.0000 18.8333 15.5500 2.42917
## 2 20.0583 2.28167 0.411250 0.161250 1.61250 25.8750 17.2083 13.7875 3.62500
## 3 23.2000 2.27833 0.559167 0.220833 1.86667 43.4167 27.2917 20.0000 8.81250
## 4 14.8125 2.35083 0.684583 0.423750 2.03750 51.5000 30.2083 26.4750 13.53330
## 5 14.2000 2.34565 0.601739 0.360000 1.76522 43.0455 26.3913 23.6696 12.93480
## 6 19.2583 2.30053 0.416250 0.235789 1.07083 25.7083 15.7083 15.2708 4.81250
## THC CO2
## 1 2.36250 444.375
## 2 2.44292 450.125
## 3 2.49917 458.417
## 4 2.77458 468.583
## 5 2.70565 464.043
## 6 2.53632 447.9173. 檢查是否有 NA 的資料。
程式碼
DATA01[!complete.cases(DATA01),]
結果
## [1] O3 CH4 CO NMHC SO2 PM10 PM2.5 NO2 NO THC CO2
## <0 列> (或零長度的 row.names)沒有缺失值問題，資料為完整，可以開始進行資料分析。
>>查看資料與繪圖
1. 基本統計量
程式碼
summary(DATA01)
結果
## O3 CH4 CO NMHC
## Min. : 3.218 Min. :1.940 Min. :0.0987 Min. :0.007083
## 1st Qu.:22.021 1st Qu.:2.044 1st Qu.:0.2346 1st Qu.:0.053333
## Median :28.654 Median :2.086 Median :0.2946 Median :0.080417
## Mean :28.800 Mean :2.113 Mean :0.3236 Mean :0.106895
## 3rd Qu.:35.283 3rd Qu.:2.172 3rd Qu.:0.3983 3rd Qu.:0.145217
## Max. :55.383 Max. :3.374 Max. :0.7512 Max. :0.423750
## SO2 PM10 PM2.5 NO2
## Min. :0.1958 Min. : 3.826 Min. : 0.913 Min. : 2.608
## 1st Qu.:0.8917 1st Qu.:13.042 1st Qu.: 6.167 1st Qu.: 5.908
## Median :1.2833 Median :20.833 Median :10.250 Median :12.508
## Mean :1.2980 Mean :23.194 Mean :12.266 Mean :12.368
## 3rd Qu.:1.6042 3rd Qu.:30.625 3rd Qu.:17.000 3rd Qu.:17.279
## Max. :3.6708 Max. :86.417 Max. :52.958 Max. :32.704
## NO THC CO2
## Min. : 0.1708 Min. :1.968 Min. :422.3
## 1st Qu.: 1.0083 1st Qu.:2.105 1st Qu.:432.0
## Median : 1.8583 Median :2.157 Median :437.6
## Mean : 3.0901 Mean :2.219 Mean :439.7
## 3rd Qu.: 3.9542 3rd Qu.:2.314 3rd Qu.:445.9
## Max. :18.7625 Max. :3.740 Max. :478.82. 變數間相關性
程式碼
cor(DATA01)
結果
## O3 CH4 CO NMHC SO2 PM10
## O3 1.00000000 -0.2653909 -0.3602427 -0.5032774 -0.1836285 -0.03374323
## CH4 -0.26539089 1.0000000 0.6299846 0.7223897 0.4227385 0.52358838
## CO -0.36024273 0.6299846 1.0000000 0.8306377 0.5256205 0.63252498
## NMHC -0.50327738 0.7223897 0.8306377 1.0000000 0.5027023 0.59177850
## SO2 -0.18362852 0.4227385 0.5256205 0.5027023 1.0000000 0.53770180
## PM10 -0.03374323 0.5235884 0.6325250 0.5917785 0.5377018 1.00000000
## PM2.5 -0.04750005 0.5756274 0.7108898 0.6252432 0.5474865 0.95764137
## NO2 -0.45871769 0.3582077 0.8172968 0.7171326 0.4916343 0.49951334
## NO -0.58227987 0.3397323 0.7197607 0.6717651 0.3300462 0.26160190
## THC -0.38509992 0.9567829 0.7605967 0.8922474 0.4874586 0.59081184
## CO2 -0.47560489 0.5614769 0.7988886 0.7567865 0.4786282 0.59764009
## PM2.5 NO2 NO THC CO2
## O3 -0.04750005 -0.4587177 -0.5822799 -0.3850999 -0.4756049
## CH4 0.57562737 0.3582077 0.3397323 0.9567829 0.5614769
## CO 0.71088980 0.8172968 0.7197607 0.7605967 0.7988886
## NMHC 0.62524317 0.7171326 0.6717651 0.8922474 0.7567865
## SO2 0.54748647 0.4916343 0.3300462 0.4874586 0.4786282
## PM10 0.95764137 0.4995133 0.2616019 0.5908118 0.5976401
## PM2.5 1.00000000 0.5331039 0.3022751 0.6388480 0.6412651
## NO2 0.53310385 1.0000000 0.7640623 0.5352192 0.6976634
## NO 0.30227512 0.7640623 1.0000000 0.5041873 0.6802579
## THC 0.63884797 0.5352192 0.5041873 1.0000000 0.6848578
## CO2 0.64126510 0.6976634 0.6802579 0.6848578 1.0000000兩兩變數間太相關需要排除!
3. 繪製散佈圖矩陣
程式碼
plot(DATA01)
結果
散佈圖矩陣（Scatterplot Matrix）是一種用於展示多個變數之間關係的圖形展示方式。它以矩陣的形式呈現，每個單元格都是一個散佈圖，顯示了兩個變數之間的關係。在散佈圖矩陣中，對角線上通常是變數的直方圖或密度圖。
散佈圖矩陣的主要目的是讓我們可以快速地視覺化多個變數之間的相互作用和相關性。
4. 變異數膨漲因子
程式碼
下載car套件才可以使用vif函數
install.packages('car') library(car)
DATA01.lm <- lm(Y ~ .,data =DATA01)
vif(DATA01.lm)
結果
## Warning: 套件 'car' 是用 R 版本 4.2.3 來建造的## 載入需要的套件：carData##
## 載入套件：'car'## 下列物件被遮斷自 'package:dplyr':
##
## recode## CH4 CO NMHC SO2 PM10 PM2.5
## 5.667093e+04 8.000079e+00 2.358244e+04 1.602005e+00 1.390121e+01 1.745729e+01
## NO2 NO THC CO2
## 4.690117e+00 3.681664e+00 1.330272e+05 3.440683e+00*判斷多元線性迴歸模型的自變數之間是否獨立，VIF值越小越好，若VIF值>10，表示自變數存在共線性，則應刪除該自變數。
所以我們應該刪除 CH4、NMHC、PM10、PM2.5、THC 這些變數(測項)
DATA01n <- DATA01[,-c(2,4,6,7,10)]
head(DATA01n)## O3 CO SO2 NO2 NO CO2
## 1 20.3625 0.384583 1.62917 15.5500 2.42917 444.375
## 2 20.0583 0.411250 1.61250 13.7875 3.62500 450.125
## 3 23.2000 0.559167 1.86667 20.0000 8.81250 458.417
## 4 14.8125 0.684583 2.03750 26.4750 13.53330 468.583
## 5 14.2000 0.601739 1.76522 23.6696 12.93480 464.043
## 6 19.2583 0.416250 1.07083 15.2708 4.81250 447.917>>建立迴歸模型
1. 建立迴歸模型
DATA01n.lm <- lm(O3 ~ CO+SO2+NO2+NO+CO2 ,data = DATA01n)2. 查看模型配適結果(含參數檢定)
程式碼
summary(DATA01n.lm)
結果
##
## Call:
## lm(formula = O3 ~ CO + SO2 + NO2 + NO + CO2, data = DATA01n)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.4604 -4.8924 -0.0238 4.0438 22.4849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 164.77882 25.42172 6.482 2.58e-10 ***
## CO 33.98468 6.43558 5.281 2.08e-07 ***
## SO2 0.26734 0.81873 0.327 0.7442
## NO2 -0.25385 0.11223 -2.262 0.0242 *
## NO -1.60102 0.19246 -8.319 1.29e-15 ***
## CO2 -0.31667 0.06066 -5.220 2.83e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.54 on 415 degrees of freedom
## Multiple R-squared: 0.3943, Adjusted R-squared: 0.387
## F-statistic: 54.03 on 5 and 415 DF, p-value: < 2.2e-16解釋: Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ’ ’ 1
是標註預測變量的顯著性”***“最顯著，然後依次減小
(由於變數顯著量過少,所以判斷為有些變數影響了其他變數的顯著性)
3. 減少變數
因此使用Forward 和 Backward(AIC準則)
Forward（前向選擇）
是一種逐步添加自變量的方法。它從一個空模型開始，然後在每一步選擇與目標變量之間具有最高相關性的一個自變量，直到滿足某個停止準則（如預先設定的顯著性水平）或沒有其他可添加的自變量為止。這種方法逐步增加自變量，每一步都根據模型的性能指標（例如，R方或AIC）來評估加入自變量後模型的改進程度。
程式碼
null = lm(Y ~ 1, data = DATA01n)
full = lm(Y ~ ., data = DATA01n)
forward.lm = step(null,scope = list(lower=null, upper=full),direction = "forward")
summary(forward.lm)
結果
## Start: AIC=1908.03
## O3 ~ 1
##
## Df Sum of Sq RSS AIC
## + NO 1 13205.8 25744 1735.7
## + CO2 1 8810.4 30139 1802.1
## + NO2 1 8195.8 30754 1810.6
## + CO 1 5054.7 33895 1851.5
## + SO2 1 1313.4 37636 1895.6
## <none> 38950 1908.0
##
## Step: AIC=1735.7
## O3 ~ NO
##
## Df Sum of Sq RSS AIC
## + CO2 1 458.26 25286 1730.1
## + CO 1 279.99 25464 1733.1
## <none> 25744 1735.7
## + NO2 1 17.87 25726 1737.4
## + SO2 1 3.20 25741 1737.7
##
## Step: AIC=1730.14
## O3 ~ NO + CO2
##
## Df Sum of Sq RSS AIC
## + CO 1 1402.90 23883 1708.1
## <none> 25286 1730.1
## + SO2 1 107.64 25178 1730.3
## + NO2 1 17.05 25268 1731.9
##
## Step: AIC=1708.11
## O3 ~ NO + CO2 + CO
##
## Df Sum of Sq RSS AIC
## + NO2 1 285.15 23597 1705.1
## <none> 23883 1708.1
## + SO2 1 0.39 23882 1710.1
##
## Step: AIC=1705.06
## O3 ~ NO + CO2 + CO + NO2
##
## Df Sum of Sq RSS AIC
## <none> 23597 1705.1
## + SO2 1 6.061 23591 1707.0##
## Call:
## lm(formula = O3 ~ NO + CO2 + CO + NO2, data = DATA01n)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.5643 -4.8498 -0.0709 4.0718 22.5923
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 163.75095 25.19897 6.498 2.33e-10 ***
## NO -1.61242 0.18906 -8.528 2.78e-16 ***
## CO2 -0.31392 0.06001 -5.231 2.67e-07 ***
## CO 34.35265 6.32934 5.428 9.73e-08 ***
## NO2 -0.24724 0.11027 -2.242 0.0255 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.532 on 416 degrees of freedom
## Multiple R-squared: 0.3942, Adjusted R-squared: 0.3883
## F-statistic: 67.66 on 4 and 416 DF, p-value: < 2.2e-16Backward（後向選擇）
是一種逐步剔除自變量的方法。它從包含所有自變量的完整模型開始，然後在每一步中，剔除對模型性能貢獻較小的一個自變量，直到滿足某個停止準則（如預先設定的顯著性水平）或沒有其他可剔除的自變量為止。這種方法逐步減少自變量，每一步都根據模型的性能指標來評估去除自變量後模型的改進程度。
程式碼
full = lm(O3 ~ ., data = DATA01n)
backward.lm = step(full, scope = list(upper=full),direction = "backward")
summary(backward.lm)
結果
## Start: AIC=1706.95
## O3 ~ CO + SO2 + NO2 + NO + CO2
##
## Df Sum of Sq RSS AIC
## - SO2 1 6.1 23597 1705.1
## <none> 23591 1707.0
## - NO2 1 290.8 23882 1710.1
## - CO2 1 1549.2 25141 1731.7
## - CO 1 1585.2 25177 1732.3
## - NO 1 3933.9 27525 1769.9
##
## Step: AIC=1705.06
## O3 ~ CO + NO2 + NO + CO2
##
## Df Sum of Sq RSS AIC
## <none> 23597 1705.1
## - NO2 1 285.1 23883 1708.1
## - CO2 1 1552.3 25150 1729.9
## - CO 1 1671.0 25268 1731.9
## - NO 1 4125.9 27723 1770.9##
## Call:
## lm(formula = O3 ~ CO + NO2 + NO + CO2, data = DATA01n)
##
## Residuals:
## Min 1Q Median 3Q Max
## -21.5643 -4.8498 -0.0709 4.0718 22.5923
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 163.75095 25.19897 6.498 2.33e-10 ***
## CO 34.35265 6.32934 5.428 9.73e-08 ***
## NO2 -0.24724 0.11027 -2.242 0.0255 *
## NO -1.61242 0.18906 -8.528 2.78e-16 ***
## CO2 -0.31392 0.06001 -5.231 2.67e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.532 on 416 degrees of freedom
## Multiple R-squared: 0.3942, Adjusted R-squared: 0.3883
## F-statistic: 67.66 on 4 and 416 DF, p-value: < 2.2e-164.精煉模型
根據原始模型、Forward 和 Backward(AIC準則的結果，可以一一推出結論。
- 先刪除不必要的
如：NO2(顯著不高移除)、NO(下面解釋)
- 新增已刪變數(是否可提高Multiple R-squared and Adjusted R-squared)
如：PM10(Multiple R-squared and Adjusted R-squared 皆有提升)
小問題1 ：這樣會顯得NO顯著不高，但刪除NO變向會使Multiple R-squared and Adjusted R-squared 降低，這樣解釋模型的能力會太低，所以還是保留NO!
小問題2 ：當NO留下則會使CO顯著性降低(>0.05)，除我們刪除CO。
最後可得出
DATA01x.lm <- lm(O3 ~ CO2+NO+PM10 ,data = DATA01)
summary(DATA01x.lm)##
## Call:
## lm(formula = O3 ~ CO2 + NO + PM10, data = DATA01)
##
## Residuals:
## Min 1Q Median 3Q Max
## -24.4916 -4.5140 0.1475 4.0131 20.8833
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 181.01648 24.94121 7.258 1.94e-12 ***
## CO2 -0.34935 0.05858 -5.963 5.28e-09 ***
## NO -1.20827 0.16030 -7.538 3.01e-13 ***
## PM10 0.22084 0.03588 6.156 1.76e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.456 on 417 degrees of freedom
## Multiple R-squared: 0.4049, Adjusted R-squared: 0.4006
## F-statistic: 94.57 on 3 and 417 DF, p-value: < 2.2e-165. 模型定型
結論
得出：O3 ~ CO2+NO+PM10
DATA01x.lm <- lm(O3 ~ CO2+NO+PM10 ,data = DATA01)
由此得知O3與 CO2、NO、PM10 這些汙染物有相關性。
圖形
# 加載所需的套件
library(ggplot2)
# 擬合回歸模型
model <- lm(O3 ~ CO2+NO+PM10 ,data = DATA01)
DATA01$predicted <- predict(model)
# 創建回歸模型圖
ggplot(DATA01, aes(x = predicted, y = O3)) +
geom_point() + # 繪制散點圖
geom_smooth(method = "lm", se = FALSE) # 添加回歸線，se = FALSE禁用置信區間## `geom_smooth()` using formula = 'y ~ x'
上面的模型中，我們使用數據集DATA01，其中CO2、NO和PM10是三個自變量，O3是因變量。然後，我們使用lm()函數擬合了一個多元線性回歸模型model，其中O3 ~ CO2+NO+PM10表示使用CO2、NO和PM10預測O3。
接下來，我們使用predict()函數計算預測值，並將預測值添加到數據集中的新列predicted。最後，使用ggplot2創建散點圖，並使用geom_smooth()函數添加擬合線。
6. 模型診斷與檢定
畫出模型診斷用的圖
程式碼
需先載入ggplot2、ggfortify才可以使用autoplot函數
library(ggplot2)
library(ggfortify)
autoplot(DATA01x.lm)
結果
## Warning: 套件 'ggfortify' 是用 R 版本 4.2.3 來建造的
Residuals vs. Fitted
表現還行的殘差圖
Normal QQ
符合常態分配
Scale-Location
這是比例位置圖，它沿 x 軸顯示回歸模型的擬合值，沿 y 軸顯示標準化殘差的平方根
然而此圖殘差的分佈在所有擬合值處大致相等
Residuals vs Leverage
這是殘差與槓桿圖，是一種診斷圖，它能夠識別回歸模型中有影響的觀察結果。
我們的Residuals vs Leverage圖中，表示我們的回歸模型中沒有任何影響點。
>>模型檢定能力
1. 殘差常態性檢定
程式碼
shapiro.test(DATA01x.lm$residual)
結果
##
## Shapiro-Wilk normality test
##
## data: DATA01x.lm$residual
## W = 0.99349, p-value = 0.06664常態性檢定的 p-value 以一般 95% 的信賴水準來說，是不拒絕虛無假設的，也就是說殘差值的常態分配假設是合理的。
2. 殘差獨立性檢定
程式碼
durbinWatsonTest(DATA01x.lm)
結果
## lag Autocorrelation D-W Statistic p-value
## 1 0.6102951 0.7757695 0
## Alternative hypothesis: rho != 0從輸出中我們可以看到檢驗統計量為 0.7757695，對應的 p 值為 0。由於此 p 值為 0，我們可以拒絕零假設並得出結論，此回歸模型中的殘差是完全正自相關的。
3. 殘差變異數同質性檢定
程式碼
ncvTest(DATA01x.lm)
結果
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.138784, Df = 1, p = 0.14362這裡發現殘差變異數同質性檢定，以一般 95% (0.05)的信賴水準來說，p = 0.14362 是不拒絕虛無假設的，也就是說殘差的變異數有符合同質性的假設。
>>結論
最後我們的迴歸模型為O3 ~ -0.34935(CO2)-1.20827(NO)+0.22084(PM10)+181.01648
由上述模型可得知當PM10下降時 O3也會隨之下降，CO2、NO上升時 O3下降，所以我們可以根據此模型來算出地區環境汙染物的平衡點，是對於該區域是最佳的!
4.4.4 子主題二
臭氧和高影響污染物的資料視覺化
4.4.5 研究主題說明
可以發現PM10的測值高時也相對造成臭氧測值的偏高，所以此主題以圖表展現，讓數據呈現減少那些汙染物可以抑制臭氧。
4.4.6 方法及結果說明
我們使用R和RMySQL套件來製作圖表，並從MySQL資料庫中擷取相關數據。以下是程式碼和結果，展示如何使用RMySQL套件連接到MySQL資料庫並製作圖表。
>>抓取資料
先從資料庫抓取所需資訊(O3&PM10)，並命名此資料表為DATA02
library(RMySQL)
# 設定資料庫連線資訊
con <- dbConnect(RMySQL::MySQL(),
user = "s10170165",
password = "s10170165@MCU1234",
host = "120.125.72.30")
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 設定工作資料庫
# 若無此資料庫，請先建立
strQuery = "USE airp15;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 執行SQL查詢並擷取數據
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL;"
DATA02 <- dbGetQuery(conn=con, statement=strQuery)
# 關閉資料庫連線
queryResult = dbDisconnect(conn=con)
#檢查資料
head(DATA02)## PM10 O3
## 1 16.8750 29.6208
## 2 16.7083 33.2708
## 3 35.2083 39.8417
## 4 20.3913 28.2750
## 5 22.6087 20.9958
## 6 15.0417 31.7667>>資料整理
這邊要以季、月做一次整理，否則在做一些圖時，以日的資料太龐大，畫出來的圖不好分析或是圖表呈現毫無意義!
抓取四季資料並建置好表。
library(RMySQL)
# 設定資料庫連線資訊
con <- dbConnect(RMySQL::MySQL(),
user = "s10170165",
password = "s10170165@MCU1234",
host = "120.125.72.30")
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 設定工作資料庫
# 若無此資料庫，請先建立
strQuery = "USE airp15;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 執行SQL查詢並擷取數據
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=3
;"
DATA_1 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=6
AND MONTH(測量日期)>=4
;"
DATA_2 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=9
AND MONTH(測量日期)>=7
;"
DATA_3 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)<=12
AND MONTH(測量日期)>=10
;"
DATA_4 <- dbGetQuery(conn=con, statement=strQuery)
# 關閉資料庫連線
queryResult = dbDisconnect(conn=con)接著整理季份資料，變成占比(季/年)
b <- colSums(DATA02)
a1 <- colSums(DATA_1)
Quarter1 <- a1 / b
a2 <- colSums(DATA_2)
Quarter2 <- a2 / b
a3 <- colSums(DATA_3)
Quarter3 <- a3 / b
a4 <- colSums(DATA_4)
Quarter4 <- a4 / b
d <- data.frame(Quarter1, Quarter2, Quarter3, Quarter4)
d## Quarter1 Quarter2 Quarter3 Quarter4
## PM10 0.2819162 0.2187419 0.2167358 0.2826061
## O3 0.2593148 0.2324895 0.2398352 0.2683605# 转置数据
e <- data.frame(t(d))
e## PM10 O3
## Quarter1 0.2819162 0.2593148
## Quarter2 0.2187419 0.2324895
## Quarter3 0.2167358 0.2398352
## Quarter4 0.2826061 0.2683605# 增加季度
Quarter <- 1:4
f <- cbind(e, Quarter)
f## PM10 O3 Quarter
## Quarter1 0.2819162 0.2593148 1
## Quarter2 0.2187419 0.2324895 2
## Quarter3 0.2167358 0.2398352 3
## Quarter4 0.2826061 0.2683605 4抓取單月份資料
library(RMySQL)
# 設定資料庫連線資訊
con <- dbConnect(RMySQL::MySQL(),
user = "s10170165",
password = "s10170165@MCU1234",
host = "120.125.72.30")
# 設定編碼
strQuery = "SET NAMES UTF8;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 設定工作資料庫
# 若無此資料庫，請先建立
strQuery = "USE airp15;"
queryResult = dbGetQuery(conn=con, statement=strQuery)
# 執行SQL查詢並擷取數據
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=1
;"
DATA_101 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=2
;"
DATA_102 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=3
;"
DATA_103 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=4
;"
DATA_104 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=5
;"
DATA_105 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=6
;"
DATA_106 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=7
;"
DATA_107 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=8
;"
DATA_108 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=9
;"
DATA_109 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=10
;"
DATA_110 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=11
;"
DATA_111 <- dbGetQuery(conn=con, statement=strQuery)
strQuery <- "SELECT `PM10`, `O3`
FROM W2022_dailyMean
WHERE 測站名稱 = ('古亭'AND'士林'AND'陽明'AND'大同'AND'松山'AND'中山'AND'萬華')
AND `O3`IS NOT NULL
AND`PM10`IS NOT NULL
AND MONTH(測量日期)=12
;"
DATA_112 <- dbGetQuery(conn=con, statement=strQuery)
# 關閉資料庫連線
queryResult = dbDisconnect(conn=con)接著整理月份資料
一月<-colSums(DATA_101)
二月<-colSums(DATA_102)
三月<-colSums(DATA_103)
四月<-colSums(DATA_104)
五月<-colSums(DATA_105)
六月<-colSums(DATA_106)
七月<-colSums(DATA_107)
八月<-colSums(DATA_108)
九月<-colSums(DATA_109)
十月<-colSums(DATA_110)
十一月<-colSums(DATA_111)
十二月<-colSums(DATA_112)
m = data.frame(一月,二月,三月,四月,五月,六月,七月,八月,九月,十月,十一月,十二月)
m## 一月 二月 三月 四月 五月 六月 七月 八月
## PM10 75556.85 53421.99 80852.89 74094.79 48673.03 40042.91 44020.20 44567.90
## O3 69644.20 66384.30 76538.31 78732.08 70773.79 41071.56 50187.53 56831.98
## 九月 十月 十一月 十二月
## PM10 72729.45 71454.73 72876.52 66013.94
## O3 89579.37 81402.93 68808.98 69769.89#轉置資料
m<-data.frame(t(m))
m## PM10 O3
## 一月 75556.85 69644.20
## 二月 53421.99 66384.30
## 三月 80852.89 76538.31
## 四月 74094.79 78732.08
## 五月 48673.03 70773.79
## 六月 40042.91 41071.56
## 七月 44020.20 50187.53
## 八月 44567.90 56831.98
## 九月 72729.45 89579.37
## 十月 71454.73 81402.93
## 十一月 72876.52 68808.98
## 十二月 66013.94 69769.89#增加季度
月<-1:12
m <- cbind(m,月)
m## PM10 O3 月
## 一月 75556.85 69644.20 1
## 二月 53421.99 66384.30 2
## 三月 80852.89 76538.31 3
## 四月 74094.79 78732.08 4
## 五月 48673.03 70773.79 5
## 六月 40042.91 41071.56 6
## 七月 44020.20 50187.53 7
## 八月 44567.90 56831.98 8
## 九月 72729.45 89579.37 9
## 十月 71454.73 81402.93 10
## 十一月 72876.52 68808.98 11
## 十二月 66013.94 69769.89 12>> 製圖
散布圖
程式碼
# 載入ggplot2套件
library(ggplot2)
# 製作散點圖
ggplot(DATA02, aes(x = PM10, y = O3)) +
geom_point() +
xlab("PM10") +
ylab("O3")結果
散點圖展示了PM10值和臭氧值之間的關係。在這個圖表中，每個散點代表一個數據點，其中x軸表示PM10值，y軸表示臭氧值。
折線圖
程式碼
# 載入ggplot2套件
library(ggplot2)
# 繪製折線圖
ggplot() +
geom_line(data=f, aes(x = Quarter, y = PM10,color = "PM10"),size=1) +
geom_point(data=f, aes(x = Quarter, y = PM10,color = "PM10"),size=3)+
geom_line(data=f, aes(x = Quarter, y = O3,color = "O3"),size=1) +
geom_point(data=f, aes(x = Quarter, y = O3,color = "O3"),size=3)+
xlab("季度")+ylab("PM10 & O3")結果
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
上述折線圖展示了PM10和臭氧值隨季節變化的趨勢。這種圖表通常用於觀察兩個變數之間的關係，並了解它們隨著時間的變化而產生的趨勢。
在折線圖中，x軸表示季度，而y軸表示數值變量，例如PM10值和臭氧值。每條折線代表一個變量的變化。如果有多個變量，它們可以使用不同的顏色或線型來區分，像這邊我們O3為紅、PM10為藍。
通過觀察折線的趨勢，我們可以得出一些初步的解釋：
趨勢方向：觀察折線的整體趨勢，我們可以發現這兩汙染物在第二、三季是處於低處，而在第一、四季處於高點。
趨勢變化率：觀察折線的斜率變化，可以了解PM10和臭氧值變化的速率。陡峭的斜率表示變化快速，而平緩的斜率表示變化緩慢。
相對位置：觀察兩條折線之間的相對位置，可以了解PM10和臭氧值之間的時間延遲關係。如果PM10的變化在臭氧值之前發生，則可以推斷PM10對臭氧的影響有一定的時間滯後。
需要注意的是，台灣PM10、O3濃度之變化除了兩者互相影響外，之變化也易受於季節影響，主要約集中於10月至2月。
柱狀圖
各別汙染物
# 載入ggplot2套件
library(ggplot2)
# 製作各別柱狀圖
ggplot(m, aes(x = 月, y = PM10)) +
geom_bar(stat = "identity", position = "dodge") +
xlab("月") +
ylab("PM10")
ggplot(m, aes(x = 月, y = O3)) +
geom_bar(stat = "identity", position = "dodge") +
xlab("月") +
ylab("O3")
在柱狀圖中，x軸表示月，而y軸個各別表示臭氧值、PM10。
可以觀察不同月份的臭氧值、PM10柱子的變化趨勢。
O3 vs PM10
# 載入ggplot2套件
library(ggplot2)
# 繪製柱狀圖
ggplot(DATA02, aes(x = PM10, y = O3)) +
geom_bar(stat = "identity", position = "dodge") +
xlab("PM10") +
ylab("O3")
可以推斷PM10對臭氧的影響是正向還是負向的。如果隨著PM10水平的增加，臭氧值的柱子也增加，則表示PM10對臭氧具有促進作用；反之，如果臭氧值的柱子隨著PM10水平的增加而減少，則表示PM10對臭氧具有抑制作用。
結論
根據研究結果，減少PM10的排放可以有效抑制臭氧（O3）的生成，這在都會區(北部)的環境管理中扮演著重要的角色。
研究選取了台北市環保署設立的七個空氣品質監測站作為研究對象，分析了2022年1至12月的空氣品質測定報告資料，並使用迴歸分析方法進行了相關性分析和探討。
研究結果顯示，在影響臭氧生成的因子中，CO2、NO和PM10對臭氧濃度的影響最大。透過迴歸分析，我們得知O3濃度的變化可以使用CO2、NO和PM10來預測，這三個因子解釋了48%的變異量（R^2 = 0.4899）。
從這些結果中，我們可以觀察到當PM10的濃度較高時，臭氧的濃度也相對較高。因此，如何有效控制這些污染物的生成量，特別是針對對臭氧生成影響最大的污染物的減排，被列為抑制臭氧生成的首要目標。
總結而言，這項研究的結果強調了北部減少PM10排放對於抑制臭氧生成的重要性。進一步的研究和環境管理措施可以針對高影響臭氧生成的污染物，尤其是PM10，以減少其排放量，從而有效改善都會區的空氣品質。
